diff --git a/examples/example09 - virtual_circuits_MASTU.ipynb b/examples/example09 - virtual_circuits_MASTU.ipynb
index d297f80..ab0918e 100644
--- a/examples/example09 - virtual_circuits_MASTU.ipynb	
+++ b/examples/example09 - virtual_circuits_MASTU.ipynb	
@@ -143,29 +143,9 @@
    "source": [
     "### Define shape targets\n",
     "\n",
-    "Next we need to define the shape targets (i.e. our quantities of interest) that we wish to monitor. There are a number of shape target pre-defined in FreeGSNKE:\n",
-    "- \"R_in\": inner midplane radius.  \n",
-    "- \"R_out\": outer midplane radius.  \n",
-    "- \"Rx_lower\": lower X-point (radial) position.\n",
-    "- \"Zx_lower\": lower X-point (vertical) position.\n",
-    "- \"Rx_upper\": upper X-point (radial) position.\n",
-    "- \"Zx_upper\": upper X-point (vertical) position.\n",
-    "- \"Rs_lower_outer\": lower strikepoint (radial) position.\n",
-    "- \"Rs_upper_outer\": upper strikepoint (radial) position.\n",
+    "Next we need to define the shape targets (i.e. our quantities of interest) that we wish to monitor. \n",
     "\n",
-    "Note that the following targets require additional options:\n",
-    "- \"Rx_lower\": approx. radial position of the lower X-point (R,Z).\n",
-    "- \"Zx_lower\": approx. vertical position of the lower X-point (R,Z).\n",
-    "- \"Rx_upper\": approx. radial position of the upper X-point (R,Z).\n",
-    "- \"Zx_upper\": approx. vertical position of the upper X-point (R,Z).\n",
-    "- \"Rs_lower_outer\": approx. (R,Z) position of the lower outer strikepoint.\n",
-    "- \"Rs_upper_outer\": approx. (R,Z) position of the upper outer strikepoint.\n",
-    "            \n",
-    "We'll see how these options are defined in a moment. While they are not all strictly required, having them is advisable as the default calculations may return spurious values in some rare cases. \n",
-    "\n",
-    "More can be added (via a feature request), though we should say that it would need to be generic enough such that its definition is well-defined across different tokamak geometries and plasmas.\n",
-    "\n",
-    "There is the option to specify **custom** shape targets if these do not work for you---we will see this shortly."
+    "The way we do this is exactly the same as was seen in Example05b with the \"plasma_descriptors\" function. The user just needs to define a function that takes in the `eq` as a (single) input and returns an array of the shape targets (sometimes called parameters) of interest. See the function defined in the next cell."
    ]
   },
   {
@@ -174,34 +154,43 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# define the targets of interest and use the VC object to calculate their values for the equilibrium above\n",
-    "targets = ['R_in', 'R_out', 'Rx_lower', 'Zx_lower', 'Rx_upper', 'Zx_upper', 'Rs_lower_outer', 'Rs_upper_outer']\n",
-    "\n",
-    "# define the target options in a dictionary (approx. (R,Z) locations of the X-points)\n",
-    "# this helps identify the correct X-point in the code\n",
-    "targets_options = dict()\n",
-    "targets_options['Rx_lower'] = np.array([0.6, -1.1])\n",
-    "targets_options['Zx_lower'] = np.array([0.6, -1.1])\n",
-    "targets_options['Rx_upper'] = np.array([0.6, 1.1])\n",
-    "targets_options['Zx_upper'] = np.array([0.6, 1.1])\n",
-    "targets_options['Rs_lower_outer'] = np.array([0.9, -1.95])\n",
-    "targets_options['Rs_upper_outer'] = np.array([0.9, 1.95])\n",
-    "\n",
-    "\n",
-    "_, target_values = VCs.calculate_targets(eq, targets, targets_options)\n",
+    "# define the descriptors function (it should return an array of values and take in an eq object)\n",
+    "def plasma_descriptors(eq):\n",
+    "\n",
+    "    # inboard/outboard midplane radii\n",
+    "    RinRout = eq.innerOuterSeparatrix()\n",
+    "\n",
+    "    # inboard midplane radius\n",
+    "    Rout = eq.innerOuterSeparatrix()[0]\n",
+    "\n",
+    "    # find lower X-point\n",
+    "    # define a \"box\" in which to search for the lower X-point\n",
+    "    XPT_BOX = [[0.33, -0.88], [0.95, -1.38]]\n",
+    "\n",
+    "    # mask those points\n",
+    "    xpt_mask = (\n",
+    "        (eq.xpt[:, 0] >= XPT_BOX[0][0])\n",
+    "        & (eq.xpt[:, 0] <= XPT_BOX[1][0])\n",
+    "        & (eq.xpt[:, 1] <= XPT_BOX[0][1])\n",
+    "        & (eq.xpt[:, 1] >= XPT_BOX[1][1])\n",
+    "    )\n",
+    "    xpts = eq.xpt[xpt_mask, 0:2].squeeze()\n",
+    "    if xpts.ndim > 1 and xpts.shape[0] > 1:\n",
+    "        opt = eq.opt[0, 0:2]\n",
+    "        dists = np.linalg.norm(xpts - opt, axis=1)\n",
+    "        idx = np.argmin(dists)  # index of closest point\n",
+    "        Rx, Zx = xpts[idx, :]\n",
+    "    else:\n",
+    "        Rx, Zx = xpts\n",
     "\n",
-    "# print\n",
-    "for i in range(len(targets)):\n",
-    "    print(targets[i] + \" = \" + str(target_values[i]))\n"
+    "    return np.array([RinRout[0], RinRout[1], Rx, Zx])"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can also define **custom** shape targets by following the template in the next cell.\n",
-    "\n",
-    "For example, here we show how to calculate the $R_{gap}$ distance in the lower divertor. This is defined (in MAST-U at least) as the radial position of the point on the separatrix at $Z=-1.5$ and the wall (on the right side). "
+    "Assign these shape targets some decriptive names (and ensure the names are in the correct order!). You can then call the function to test it works correctly. "
    ]
   },
   {
@@ -210,67 +199,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# any new target function should take the equilibrium object as input\n",
-    "import shapely as sh # requires the shapely package\n",
-    "\n",
-    "def R_gap_lower(eq):\n",
-    "    \n",
-    "    # find contour object for psi_boundary\n",
-    "    if eq._profiles.flag_limiter:\n",
-    "        cs = plt.contour(\n",
-    "            eq.R, eq.Z, eq.psi(), levels=[eq._profiles.psi_bndry]\n",
-    "        )\n",
-    "    else:\n",
-    "        cs = plt.contour(\n",
-    "            eq.R, eq.Z, eq.psi(), levels=[eq._profiles.xpt[0][2]]\n",
-    "        )\n",
-    "    plt.close()  # this isn't the most elegant but we don't need the plot iteq\n",
-    "\n",
-    "    # for each item in the contour object there's a list of points in (r,z) (i.e. a line)\n",
-    "    psi_boundary_lines = []\n",
-    "    for i, item in enumerate(cs.allsegs[0]):\n",
-    "        psi_boundary_lines.append(item)\n",
-    "\n",
-    "    # use the shapely package to find where each psi_boundary_line intersects the z line\n",
-    "    gaps = []\n",
-    "    z_line = np.array([[0.5,-1.5],[0.8,-1.5]])\n",
-    "    curve1 = sh.LineString(z_line)\n",
-    "    for j, line in enumerate(psi_boundary_lines):\n",
-    "        curve2 = sh.LineString(line)\n",
-    "\n",
-    "        # find the intersection point(s)\n",
-    "        intersection = curve2.intersection(curve1)\n",
-    "\n",
-    "        # extract intersection point(s)\n",
-    "        if intersection.geom_type == \"Point\":\n",
-    "            gaps.append(np.squeeze(np.array(intersection.xy).T))\n",
-    "        elif intersection.geom_type == \"MultiPoint\":\n",
-    "            gaps.append(\n",
-    "                np.squeeze(\n",
-    "                    np.array([geom.xy for geom in intersection.geoms])\n",
-    "                )\n",
-    "            )\n",
-    "\n",
-    "    gap_point = np.array(gaps).squeeze()\n",
-    "\n",
-    "    return gap_point[0] # select R position\n",
-    "\n",
-    "# we then create a list of lists to store the target name and its function, e.g. [[\"new_target_name\", ...], [new_target_function(eq),...]]\n",
-    "non_standard_targets = [[\"R_gap_lower\"], [R_gap_lower]]\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# we can now include our custom target functions\n",
-    "all_target_names, all_target_values = VCs.calculate_targets(eq, targets, targets_options, non_standard_targets)\n",
-    "\n",
-    "# print\n",
-    "for i in range(len(all_target_names)):\n",
-    "    print(all_target_names[i] + \" = \" + str(all_target_values[i]))\n"
+    "target_names = [\"Rin\", \"Rout\", \"Rx\", \"Zx\"]\n",
+    "target_values = plasma_descriptors(eq)"
    ]
   },
   {
@@ -293,17 +223,13 @@
     "eq.tokamak.plot(axis=ax1, show=False)\n",
     "ax1.plot(eq.tokamak.wall.R, eq.tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n",
     "\n",
-    "ax1.scatter(all_target_values[0], 0.0, s=100, color='green', marker='x', zorder=20, label=all_target_names[0])\n",
-    "ax1.scatter(all_target_values[1], 0.0, s=100, color='blue', marker='x', zorder=20, label=all_target_names[1])\n",
-    "ax1.scatter(all_target_values[2], all_target_values[3], s=100, color='m', marker='*', zorder=20, label=f\"({all_target_names[2]},{all_target_names[3]})\")\n",
-    "ax1.scatter(all_target_values[4], all_target_values[5], s=100, color='k', marker='*', zorder=20, label=f\"({all_target_names[4]},{all_target_names[5]})\")\n",
-    "ax1.scatter(all_target_values[6], -1.95 , s=100, color='k', marker='x', zorder=20, label=f\"{all_target_names[6]}\")\n",
-    "ax1.scatter(all_target_values[7], 1.95 , s=100, color='r', marker='x', zorder=20, label=f\"{all_target_names[7]}\")\n",
-    "ax1.scatter(all_target_values[8], -1.5 , s=100, color='orange', marker='o', zorder=5, label=f\"{all_target_names[8]}\")\n",
+    "ax1.scatter(target_values[0], 0.0, s=100, color='green', marker='x', zorder=20, label=target_names[0])\n",
+    "ax1.scatter(target_values[1], 0.0, s=100, color='blue', marker='x', zorder=20, label=target_names[1])\n",
+    "ax1.scatter(target_values[2], target_values[3], s=100, color='m', marker='*', zorder=20, label=f\"({target_names[2]},{target_names[3]})\")\n",
     "\n",
     "ax1.set_xlim(0.1, 2.15)\n",
     "ax1.set_ylim(-2.25, 2.25)\n",
-    "ax1.legend(loc=\"center\")\n",
+    "ax1.legend(loc=\"upper right\")\n",
     "plt.tight_layout()"
    ]
   },
@@ -374,21 +300,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# define which coils we wish to calculate the shape (finite difference) derivatives (and therefore VCs) for\n",
-    "coils = eq.tokamak.coils_list[0:12]\n",
-    "print(coils)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "# next we need to define which coils we wish to calculate the shape (finite difference) derivatives (and therefore VCs) for\n",
+    "print(\"All active coils :\", eq.tokamak.coils_list[0:12])\n",
+    "\n",
     "# here we'll look at a subset of the coils and the targets\n",
-    "coils = ['D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7', 'P4', 'P5']\n",
-    "targets = ['R_in', 'R_out', 'Rx_lower', 'Zx_lower', 'Rs_lower_outer']\n",
-    "non_standard_targets = [[\"R_gap_lower\"], [R_gap_lower]]"
+    "coils = ['PX', 'D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7', 'P4', 'P5']\n",
+    "print(\"Coils to calculate VCs for :\", coils)"
    ]
   },
   {
@@ -397,19 +314,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# calculate the shape and VC matrices as follows\n",
+    "# calculate the shape and VC matrices using the in-built method\n",
     "VCs.calculate_VC(\n",
     "    eq=eq,\n",
     "    profiles=profiles,\n",
     "    coils=coils,\n",
-    "    targets=targets,\n",
-    "    targets_options=targets_options,\n",
-    "    target_dIy=1e-3,\n",
-    "    non_standard_targets=non_standard_targets,\n",
-    "    starting_dI=None,\n",
+    "    target_names=target_names,\n",
+    "    target_calculator=plasma_descriptors,\n",
+    "    starting_dI=None,    \n",
     "    min_starting_dI=50,\n",
     "    verbose=True,\n",
-    "    VC_name=\"VC_for_lower_targets\", # name for the VC\n",
+    "    name=\"VC_for_lower_targets\", # name for the VC\n",
     "    )"
    ]
   },
@@ -446,9 +361,7 @@
    "outputs": [],
    "source": [
     "# let's remind ourselves of the targets\n",
-    "print(targets)\n",
-    "for i in range(len(non_standard_targets[0])):\n",
-    "    print(non_standard_targets[0][i])"
+    "print(target_names)"
    ]
   },
   {
@@ -458,7 +371,7 @@
    "outputs": [],
    "source": [
     "# let's set the requested shifts (units are metres) for the full set of targets\n",
-    "all_requested_target_shifts = [0.01, -0.01, 0.0, 0.0, 0.0, 0.0]"
+    "requested_target_shifts = [0.05, -0.05, 0.0, 0.0]"
    ]
   },
   {
@@ -479,11 +392,11 @@
    "outputs": [],
    "source": [
     "# we apply the VCs to the desired shifts, using the VC_object we just created\n",
-    "eq_new, profiles_new, all_target_names, new_target_values, old_target_values = VCs.apply_VC(\n",
+    "eq_new, profiles_new, new_target_values, old_target_values = VCs.apply_VC(\n",
     "    eq=eq,\n",
     "    profiles=profiles,\n",
     "    VC_object=VCs.VC_for_lower_targets,\n",
-    "    all_requested_target_shifts=all_requested_target_shifts,\n",
+    "    requested_target_shifts=requested_target_shifts,\n",
     "    verbose=True,\n",
     ")"
    ]
@@ -501,7 +414,7 @@
     "eq.tokamak.plot(axis=ax1, show=False)\n",
     "\n",
     "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n",
-    "ax1.contour(eq.R, eq.Z, eq_new.psi(), levels=[eq_new.psi_bndry], colors='b', linestyles=\"--\")\n",
+    "ax1.contour(eq.R, eq.Z, eq_new.psi(), levels=[eq_new.psi_bndry], colors='b')\n",
     "\n",
     "ax1.set_xlim(0.1, 2.15)\n",
     "ax1.set_ylim(-2.25, 2.25)\n",
@@ -523,8 +436,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "for i in range(len(all_target_names)):\n",
-    "    print(f\"Difference in {all_target_names[i]} = {np.round(new_target_values[i] - old_target_values[i],3)} vs. requested = {(all_requested_target_shifts)[i]}.\")\n"
+    "for i in range(len(target_names)):\n",
+    "    print(f\"Difference in {target_names[i]} = {np.round(new_target_values[i] - old_target_values[i],3)} vs. requested = {(requested_target_shifts)[i]}.\")"
    ]
   },
   {
@@ -546,15 +459,15 @@
     "fig1, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(12,12), dpi=80)\n",
     "\n",
     "ax1.grid(True, which='both', alpha=0.5)\n",
-    "ax1.scatter(all_target_names, all_requested_target_shifts, color='red', marker='o', s=150, label=\"Requested\")\n",
-    "ax1.scatter(all_target_names, np.array(new_target_values) - np.array(old_target_values), color='royalblue', marker='o', s=75, label=\"Actual\")\n",
+    "ax1.scatter(target_names, requested_target_shifts, color='red', marker='o', s=150, label=\"Requested\")\n",
+    "ax1.scatter(target_names, np.array(new_target_values) - np.array(old_target_values), color='royalblue', marker='o', s=75, label=\"Actual\")\n",
     "ax1.set_xlabel(\"Target\")\n",
     "ax1.set_ylabel(\"Shift [m]\")\n",
     "ax1.legend()\n",
     "# ax1.set_ylim([-max(targets_shift+non_standard_targets_shift)*1.1, max(targets_shift+non_standard_targets_shift)*1.1])\n",
     "\n",
     "ax2.grid(True, which='both', alpha=0.5)\n",
-    "ax2.scatter(all_target_names, rel_diff, color='red', marker='o', s=150, edgecolors='black', label=\"Requested\")\n",
+    "ax2.scatter(target_names, rel_diff, color='red', marker='o', s=150, edgecolors='black', label=\"Requested\")\n",
     "ax2.set_xlabel(\"Target\")\n",
     "ax2.set_ylabel(\"Relative shift\")\n",
     "labels = ax2.get_xticklabels()\n",
@@ -570,9 +483,9 @@
    "source": [
     "### Alternative VCs\n",
     "\n",
-    "Note that when the `VC_calculate` method is called, a `VC_name` can be provided. This will store the key data in a subclass with `VC_name`. This useful for recalling from which targets and coils a given shape matrix or VC matrix was calculated. \n",
+    "Note that when the `VC_calculate` method is called, a `name` can be provided. This will store the key data in a subclass with `name`. This useful for recalling from which targets and coils a given shape matrix or VC matrix was calculated. \n",
     "\n",
-    "Note that if no `VC_name` is provided the default is used (\"latest_VC\")."
+    "Note that if no `name` is provided the default is used (\"latest_VC\")."
    ]
   },
   {
@@ -598,28 +511,57 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# choose new coils and targets\n",
-    "coils = ['D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7', 'P4', 'P5']\n",
-    "targets = ['Rx_upper', 'Zx_upper', 'Rs_upper_outer'] # this time we use the upper targets\n",
-    "targets_options = dict()\n",
-    "targets_options['Rx_upper'] = np.array([0.6, 1.1])\n",
-    "targets_options['Zx_upper'] = np.array([0.6, 1.1])\n",
-    "targets_options['Rs_upper_outer'] = np.array([1.15, 2.1])\n",
+    "# define the descriptors function (it should return an array of values and take in an eq object)\n",
+    "def plasma_descriptors_new(eq):\n",
+    "\n",
+    "    # inboard/outboard midplane radii\n",
+    "    RinRout = eq.innerOuterSeparatrix()\n",
+    "\n",
+    "    # find lower X-point\n",
+    "    # define a \"box\" in which to search for the lower X-point\n",
+    "    XPT_BOX = [[0.33, -0.88], [0.95, -1.38]]\n",
+    "\n",
+    "    # mask those points\n",
+    "    xpt_mask = (\n",
+    "        (eq.xpt[:, 0] >= XPT_BOX[0][0])\n",
+    "        & (eq.xpt[:, 0] <= XPT_BOX[1][0])\n",
+    "        & (eq.xpt[:, 1] <= XPT_BOX[0][1])\n",
+    "        & (eq.xpt[:, 1] >= XPT_BOX[1][1])\n",
+    "    )\n",
+    "    xpts = eq.xpt[xpt_mask, 0:2].squeeze()\n",
+    "    if xpts.ndim > 1 and xpts.shape[0] > 1:\n",
+    "        opt = eq.opt[0, 0:2]\n",
+    "        dists = np.linalg.norm(xpts - opt, axis=1)\n",
+    "        idx = np.argmin(dists)  # index of closest point\n",
+    "        Rx, Zx = xpts[idx, :]\n",
+    "    else:\n",
+    "        Rx, Zx = xpts\n",
+    "\n",
+    "    return np.array([RinRout[0], Rx, Zx])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# this time we use only the divertor coils and a subset of targets\n",
+    "coils = [\"PX\", 'D1', 'D2', 'D3', 'Dp', 'D5', 'D6', 'D7']\n",
+    "target_names =['Rin', 'Rx', 'Zx']\n",
     "\n",
     "\n",
-    "# calculate the shape and VC matrices\n",
+    "# calculate the shape and VC matrices as follows\n",
     "VCs.calculate_VC(\n",
     "    eq=eq,\n",
     "    profiles=profiles,\n",
     "    coils=coils,\n",
-    "    targets=targets,\n",
-    "    targets_options=targets_options,\n",
-    "    target_dIy=1e-3,\n",
-    "    non_standard_targets=None,\n",
+    "    target_names=target_names,\n",
+    "    target_calculator=plasma_descriptors_new,\n",
     "    starting_dI=None,\n",
     "    min_starting_dI=50,\n",
     "    verbose=True,\n",
-    "    VC_name=\"VC_for_upper_targets\",\n",
+    "    name=\"VC_for_upper_targets\", # name for the VC\n",
     "    )"
    ]
   },
@@ -647,27 +589,18 @@
    "outputs": [],
    "source": [
     "# we'll shift the upper strikepoint (R value) \n",
-    "all_requested_target_shifts = [0.0, 0.0, 0.02]\n",
+    "requested_target_shifts = [0.0, -0.025, 0.03]\n",
     "\n",
     "# we apply the VCs to the some desired shifts, using the VC_object we just created\n",
-    "eq_new, profiles_new, target_names, new_target_values, old_target_values = VCs.apply_VC(\n",
+    "eq_new, profiles_new, new_target_values, old_target_values = VCs.apply_VC(\n",
     "    eq=eq,\n",
     "    profiles=profiles,\n",
     "    VC_object=VCs.VC_for_upper_targets,\n",
-    "    all_requested_target_shifts=all_requested_target_shifts,\n",
+    "    requested_target_shifts=requested_target_shifts,\n",
     "    verbose=True,\n",
     ")"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "VCs.VC_for_upper_targets.targets"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -681,7 +614,7 @@
     "eq.tokamak.plot(axis=ax1, show=False)\n",
     "\n",
     "ax1.contour(eq.R, eq.Z, eq.psi(), levels=[eq.psi_bndry], colors='r')\n",
-    "ax1.contour(eq.R, eq.Z, eq_new.psi(), levels=[eq_new.psi_bndry], colors='b', linestyles=\"--\")\n",
+    "ax1.contour(eq.R, eq.Z, eq_new.psi(), levels=[eq_new.psi_bndry], colors='b')\n",
     "\n",
     "ax1.set_xlim(0.1, 2.15)\n",
     "ax1.set_ylim(-2.25, 2.25)\n",
@@ -701,7 +634,7 @@
     "fig1, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(12,12), dpi=80)\n",
     "\n",
     "ax1.grid(True, which='both', alpha=0.5)\n",
-    "ax1.scatter(target_names, all_requested_target_shifts, color='red', marker='o', s=150, label=\"Requested\")\n",
+    "ax1.scatter(target_names, requested_target_shifts, color='red', marker='o', s=150, label=\"Requested\")\n",
     "ax1.scatter(target_names, np.array(new_target_values) - np.array(old_target_values), color='royalblue', marker='o', s=75, label=\"Actual\")\n",
     "ax1.set_xlabel(\"Target\")\n",
     "ax1.set_ylabel(\"Shift [m]\")\n",
diff --git a/freegsnke/virtual_circuits.py b/freegsnke/virtual_circuits.py
index be0d1d3..73358b9 100644
--- a/freegsnke/virtual_circuits.py
+++ b/freegsnke/virtual_circuits.py
@@ -19,7 +19,7 @@
 along with FreeGSNKE.  If not, see <http://www.gnu.org/licenses/>. 
 """
 
-from copy import deepcopy
+from datetime import datetime
 
 import numpy as np
 
@@ -37,11 +37,9 @@ def __init__(
         profiles,
         shape_matrix,
         VCs_matrix,
-        targets,
-        targets_val,
-        targets_options,
-        non_standard_targets,
+        target_names,
         coils,
+        target_calculator,
     ):
         """
         Store the key quantities from the VirtualCircuitHandling calculations.
@@ -60,19 +58,12 @@ def __init__(
         VCs_matrix : np.array
             The array storing the VCs between the targets and coils given in 'targets'
             and 'coils'.
-        targets : list
-            The list of targets used to calculate the shape_matrix and VCs_matrix.
-        targets_val : np.array
-            The array of target values.
-        targets_options : dict
-            Dictionary of additional parameters required to calculate the
-            'targets'.
-        non_standard_targets : list
-            List of lists of additional (non-standard) target functions to use. Takes the
-            form [["new_target_name",...], [function(eq),...]], where function calcualtes the target
-            using the eq object.
+        target_names : list
+            A list of the target names, e.g [Rin, Rout, Rx, Zx, ...] (must be same length as array from target_calculator).
         coils : list
             The list of coils used to calculate the shape_matrix and VCs_matrix.
+        target_calculator : function
+            Function returning an array of the shape targets (VC will be calculated for ALL of these targets).
         """
 
         self.name = name
@@ -80,15 +71,9 @@ def __init__(
         self.profiles = profiles
         self.shape_matrix = shape_matrix
         self.VCs_matrix = VCs_matrix
-        self.targets = targets
-        self.targets_val = targets_val
-        self.non_standard_targets = non_standard_targets
+        self.target_names = target_names
         self.coils = coils
-        self.targets_options = targets_options
-        if self.non_standard_targets is not None:
-            self.len_non_standard_targets = len(self.non_standard_targets[0])
-        else:
-            self.len_non_standard_targets = 0
+        self.target_calculator = target_calculator
 
 
 class VirtualCircuitHandling:
@@ -108,7 +93,7 @@ def __init__(
         """
 
         # name to store the VC under
-        self.default_VC_name = "latest_VC"
+        self.default_VC_name = f"VC_{datetime.today().strftime('%Y%m%d')}"
 
     def define_solver(self, solver, target_relative_tolerance=1e-7):
         """
@@ -130,172 +115,6 @@ def define_solver(self, solver, target_relative_tolerance=1e-7):
         self.solver = solver
         self.target_relative_tolerance = target_relative_tolerance
 
-    def calculate_targets(
-        self, eq, targets, targets_options=None, non_standard_targets=None
-    ):
-        """
-        For the given equilibrium, this function calculates the targets
-        specified in the targets list.
-
-        Parameters
-        ----------
-        eq : object
-            The equilibrium object.
-        targets : list
-            List of strings containing the targets of interest. Currently supported targets
-            are:
-            - "R_in": inner midplane radius.
-            - "R_out": outer midplane radius.
-            - "Rx_lower": lower X-point (radial) position.
-            - "Zx_lower": lower X-point (vertical) position.
-            - "Rx_upper": upper X-point (radial) position.
-            - "Zx_upper": upper X-point (vertical) position.
-            - "Rs_lower_outer": lower strikepoint (radial) position.
-            - "Rs_upper_outer": upper strikepoint (radial) position.
-        targets_options : dict
-            Dictionary of additional parameters required to calculate the
-            'targets'. Options are required for:
-            - "Rx_lower": approx. (R,Z) position of the lower X-point.
-            - "Zx_lower": approx. (R,Z) position of the lower X-point.
-            - "Rx_upper": approx. (R,Z) position of the upper X-point.
-            - "Zx_upper": approx. (R,Z) position of the upper X-point.
-            - "Rs_lower_outer": approx. (R,Z) position of the lower outer strikepoint.
-            - "Rs_upper_outer": approx. (R,Z) position of the upper outer strikepoint.
-        non_standard_targets : list
-            List of lists of additional (non-standard) target functions to use. Takes the
-            form [["new_target_name",...], [function(eq),...]], where function calcualtes the target
-            using the eq object.
-
-        Returns
-        -------
-        list
-            Returns the original list of targets (plus any additional tagets specified
-            in non_standard_targets).
-        np.array
-            Returns a 1D array of the target values (in same order as 'targets' input).
-        """
-
-        def get_xpoint_coordinate(target_name, coord_index, use_max_z=False):
-            """
-            Extract X-point coordinate (R or Z).
-
-            Parameters
-            ----------
-            target_name : str
-                Name of target (for error messages and options lookup)
-            coord_index : int
-                0 for R, 1 for Z
-            use_max_z : bool
-                If True, use argmax for fallback (upper), else argmin (lower)
-            """
-            if target_name in targets_options:
-                loc = targets_options[target_name]
-                xpts = eq.xpt[:, 0:2]
-                x_point_ind = np.argmin(np.sum((xpts - loc) ** 2, axis=1))
-                return xpts[x_point_ind, coord_index]
-            else:
-                print(
-                    f"Use of the 'target_options' input for {target_name} is advised!"
-                )
-                xpts = eq.xpt[0:2, 0:2]
-                x_point_ind = (
-                    np.argmax(xpts[:, 1]) if use_max_z else np.argmin(xpts[:, 1])
-                )
-                return xpts[x_point_ind, coord_index]
-
-        def get_strikepoint_r(target_name, upper=False):
-            """Extract strikepoint R coordinate."""
-            if target_name in targets_options:
-                loc = targets_options[target_name]
-                strikes = eq.strikepoints()
-                strike_ind = np.argmin(np.sum((strikes - loc) ** 2, axis=1))
-                return strikes[strike_ind, 0]
-            else:
-                print(
-                    f"Use of the 'target_options' input for {target_name} is advised!"
-                )
-                strikes = eq.strikepoints()
-                if strikes.shape[0] > 4:
-                    print(
-                        f"More than four strikepoints located, use of 'target_options' "
-                        f"input for {target_name} is strongly advised!"
-                    )
-                # Filter by Z sign and find max R
-                filtered = (
-                    strikes[strikes[:, 1] > 0] if upper else strikes[strikes[:, 1] < 0]
-                )
-                return filtered[np.argmax(filtered[:, 0]), 0]
-
-        # Initialise
-        rinout_flag = False  # turns True when calculated
-        if targets_options is None:
-            targets_options = {}
-
-        # output targets
-        final_targets = deepcopy(targets)
-        n_targets = len(targets) + (
-            len(non_standard_targets[0]) if non_standard_targets else 0
-        )
-        target_vec = np.zeros(n_targets)
-
-        # Handle non-standard targets only case
-        if len(targets) == 0 and non_standard_targets is not None:
-            final_targets += non_standard_targets[0]
-            for j, func in enumerate(non_standard_targets[1]):
-                target_vec[j] = func(eq)
-            return final_targets, target_vec
-
-        # Calculate standard targets
-        for i, target in enumerate(targets):
-            if target == "R_in":
-                if not rinout_flag:
-                    rin, rout = eq.innerOuterSeparatrix()
-                    rinout_flag = True
-                target_vec[i] = rin
-
-            elif target == "R_out":
-                if not rinout_flag:
-                    rin, rout = eq.innerOuterSeparatrix()
-                    rinout_flag = True
-                target_vec[i] = rout
-
-            elif target == "Rx_lower":
-                target_vec[i] = get_xpoint_coordinate(
-                    target, coord_index=0, use_max_z=False
-                )
-
-            elif target == "Zx_lower":
-                target_vec[i] = get_xpoint_coordinate(
-                    target, coord_index=1, use_max_z=False
-                )
-
-            elif target == "Rx_upper":
-                target_vec[i] = get_xpoint_coordinate(
-                    target, coord_index=0, use_max_z=True
-                )
-
-            elif target == "Zx_upper":
-                target_vec[i] = get_xpoint_coordinate(
-                    target, coord_index=1, use_max_z=True
-                )
-
-            elif target == "Rs_lower_outer":
-                target_vec[i] = get_strikepoint_r(target, upper=False)
-
-            elif target == "Rs_upper_outer":
-                target_vec[i] = get_strikepoint_r(target, upper=True)
-
-            else:
-                raise ValueError(f"Undefined target: {target}.")
-
-        # Add non-standard targets
-        if non_standard_targets is not None:
-            final_targets += non_standard_targets[0]
-            for j, func in enumerate(non_standard_targets[1]):
-                target_vec[len(targets) + j] = func(eq)
-
-        return final_targets, target_vec
-
     def build_current_vec(self, eq, coils):
         """
         For the given equilibrium, this function stores the coil currents
@@ -370,11 +189,15 @@ def assign_currents_solve_GS(self, currents_vec, coils, target_relative_toleranc
         self.assign_currents(currents_vec, coils, eq=self._eq2)
 
         # solve for equilibrium
-        self.solver.forward_solve(
-            self._eq2,
-            self._profiles2,
-            target_relative_tolerance=target_relative_tolerance,
-        )
+        try:
+            self.solver.forward_solve(
+                self._eq2,
+                self._profiles2,
+                target_relative_tolerance=target_relative_tolerance,
+                # suppress=True,
+            )
+        except AttributeError:
+            raise AttributeError("Solver not defined. Call define_solver() first.")
 
     def prepare_build_dIydI_j(
         self, j, coils, target_dIy, starting_dI, min_curr=1e-4, max_curr=300
@@ -422,7 +245,13 @@ def prepare_build_dIydI_j(
         dIy_0 = self._eq2.limiter_handler.Iy_from_jtor(self._profiles2.jtor) - self.Iy
 
         # relative norm of plasma current change
-        rel_ndIy_0 = np.linalg.norm(dIy_0) / self._nIy
+        norm_dIy_0 = np.linalg.norm(dIy_0)
+        if norm_dIy_0 < 1e-10:
+            raise ZeroDivisionError(
+                "Norm of change in jtor is near-zero for this Jacobian, please increase 'starting_dI' parameter."
+            )
+        else:
+            rel_ndIy_0 = norm_dIy_0 / self._nIy
 
         # scale the starting_dI to match the target
         final_dI = starting_dI * target_dIy / rel_ndIy_0
@@ -437,9 +266,6 @@ def build_dIydI_j(
         self,
         j,
         coils,
-        targets,
-        targets_options,
-        non_standard_targets=None,
         verbose=False,
     ):
         """
@@ -456,14 +282,6 @@ def build_dIydI_j(
             Index identifying the current to be varied. Indexes as in self.currents_vec.
         coils : list
             List of strings containing the names of the coil currents to be assigned.
-        targets : list
-            List of strings containing the targets of interest. See above for supported targets.
-        targets_options : dict
-            Dictionary of additional parameters required to calculate the 'targets' (see above).
-        non_standard_targets : list
-            List of lists of additional (non-standard) target functions to use. Takes the
-            form [["new_target_name",...], [function(eq),...]], where function calcualtes the target
-            using the eq object.
         verbose: bool
             Display output (or not).
 
@@ -473,6 +291,10 @@ def build_dIydI_j(
             VC object modifed in place.
         """
 
+        # print some output
+        if verbose:
+            print(f"Coil {coils[j]}")
+
         # store dI
         final_dI = 1.0 * self.final_dI_record[j]
 
@@ -486,23 +308,9 @@ def build_dIydI_j(
         self.assign_currents_solve_GS(currents, coils, self.target_relative_tolerance)
 
         # calculate finite difference of targets wrt to the coil current
-        _, self._target_vec_1 = self.calculate_targets(
-            self._eq2, targets, targets_options, non_standard_targets
-        )
-        dtargets = self._target_vec_1 - self._targets_vec
-        # self._dtargetsdIj = dtargets / final_dI
+        self._target_vec_1 = self.target_calculator(self._eq2)
 
-        # print some output
-        if verbose:
-            print(f"{j}th coil ({coils[j]}) using scaled current shift {final_dI}.")
-            # print(
-            #     "Direction (coil)",
-            #     j,
-            #     ", gradient calculated on the finite difference: norm(deltaI) = ",
-            #     final_dI,
-            #     ", norm(deltaIy) =",
-            #     np.linalg.norm(dIy_1),
-            # )
+        dtargets = self._target_vec_1 - self._targets_vec
 
         return dtargets / final_dI
 
@@ -511,14 +319,13 @@ def calculate_VC(
         eq,
         profiles,
         coils,
-        targets,
-        targets_options,
-        non_standard_targets=None,
+        target_names,
+        target_calculator,
         target_dIy=1e-3,
         starting_dI=None,
         min_starting_dI=50,
         verbose=False,
-        VC_name=None,
+        name=None,
     ):
         """
         Calculate the "virtual circuits" matrix:
@@ -540,23 +347,19 @@ def calculate_VC(
             The profiles object.
         coils : list
             List of strings containing the names of the coil currents to be assigned.
-        targets : list
-            List of strings containing the targets of interest. See above for supported targets.
-        targets_options : dict
-            Dictionary of additional parameters required to calculate the 'targets' (see above).
-        non_standard_targets : list
-            List of lists of additional (non-standard) target functions to use. Takes the
-            form [["new_target_name",...], [function(eq),...]], where function calcualtes the target
-            using the eq object.
+        target_names : list
+            A list of the target names, e.g [Rin, Rout, Rx, Zx, ...] (must be same length as array from target_calculator).
+        target_calculator : function
+            Function returning an array of the shape targets (VC will be calculated for ALL of these targets).
         target_dIy : float
             Target value for the norm of delta(I_y), from which the finite difference derivative is calculated.
-        starting_dI : float
-            Initial value to be used as delta(I_j) to infer the slope of norm(delta(I_y))/delta(I_j).
+        starting_dI : array
+            Initial current perturbations [Amps] to be used as delta(I_j) to infer the slope of norm(delta(I_y))/delta(I_j).
         min_starting_dI : float
             Minimum starting_dI value to be used as delta(I_j): to infer the slope of norm(delta(I_y))/delta(I_j).
         verbose: bool
             Display output (or not).
-        VC_name: str
+        name: str
             Name to store the VC under (in the 'VirtualCircuit' class).
 
         Returns
@@ -569,22 +372,37 @@ def calculate_VC(
         # store original currents
         self.build_current_vec(eq, coils)
 
+        # store function to calculate targets from equilibrium
+        if target_calculator is None:
+            raise ValueError("You need to input a 'target_calculator' function!")
+        self.target_calculator = target_calculator
+
+        # calculate the targets from the equilibrium
+        self._targets_vec = self.target_calculator(eq)
+
+        if target_names is None:
+            raise ValueError("You need to input a list of 'target_names'!")
+        elif len(target_names) != len(self._targets_vec):
+            raise ValueError(
+                "Number of 'target_names' does not match length of array from 'target_calculator' function!"
+            )
+        self.target_names = target_names
+
         # solve static GS problem (it's already solved?)
-        self.solver.forward_solve(
-            eq=eq,
-            profiles=profiles,
-            target_relative_tolerance=self.target_relative_tolerance,
-        )
+        try:
+            self.solver.forward_solve(
+                eq=eq,
+                profiles=profiles,
+                target_relative_tolerance=self.target_relative_tolerance,
+                suppress=True,
+            )
+        except AttributeError:
+            raise AttributeError("Solver not defined. Call define_solver() first.")
 
         # store the flattened plasma current vector (and its norm)
         self.Iy = eq.limiter_handler.Iy_from_jtor(profiles.jtor).copy()
         self._nIy = np.linalg.norm(self.Iy)
 
-        # calculate the targets from the equilibrium
-        targets_new, self._targets_vec = self.calculate_targets(
-            eq, targets, targets_options, non_standard_targets
-        )
-
         # define starting_dI using currents if not given
         if starting_dI is None:
             starting_dI = np.abs(self.currents_vec.copy()) * target_dIy
@@ -593,11 +411,13 @@ def calculate_VC(
             )
 
         if verbose:
-            print("---")
-            print("Preparing the scaled current shifts with respect to the:")
+            print("--- Stage one ---")
+            print(
+                f"Re-sizing each initial coil current shift so that it produces a {np.round(target_dIy*100,2)}% change in plasma current density from the input equilibrium."
+            )
 
         # storage matrices
-        shape_matrix = np.zeros((len(targets_new), len(coils)))
+        shape_matrix = np.zeros((len(self._targets_vec), len(coils)))
         self.final_dI_record = np.zeros(len(coils))
 
         # make copies of the newly solved equilibrium and profile objects
@@ -608,55 +428,53 @@ def calculate_VC(
         # for each coil, prepare by inferring delta(I_j) corresponding to a change delta(I_y)
         # with norm(delta(I_y)) = target_dIy
         for j in np.arange(len(coils)):
+            self.prepare_build_dIydI_j(j, coils, target_dIy, starting_dI[j])
             if verbose:
                 print(
-                    f"{j}th coil ({coils[j]}) using initial current shift {starting_dI[j]}."
+                    f"Coil {coils[j]} (original current shift = {np.round(starting_dI[j],2)} [A] --> scaled current shift {np.round(self.final_dI_record[j],2)} [A])."
                 )
-            self.prepare_build_dIydI_j(j, coils, target_dIy, starting_dI[j])
 
         if verbose:
-            print("---")
-            print("Building the shape matrix with respect to the:")
+            print("--- Stage two ---")
+            print(
+                "Building the shape matrix (Jacobian) of the shape parameter changes wrt scaled current shifts for each coil:"
+            )
 
         # for each coil, build the Jacobian using the value of delta(I_j) inferred earlier
         # by self.prepare_build_dIydI_j.
         for j in np.arange(len(coils)):
             # each shape matrix row is derivative of targets wrt the final coil current change
-            shape_matrix[:, j] = self.build_dIydI_j(
-                j, coils, targets, targets_options, non_standard_targets, verbose
-            )
+            shape_matrix[:, j] = self.build_dIydI_j(j, coils, verbose)
 
         # store the data in its own (new) class
-        if VC_name is None:
-            VC_name = self.default_VC_name
+        if name is None:
+            name = self.default_VC_name
+
+        print("--- Stage three ---")
+        print("Inverting the shape matrix to get the virtual circuit matrix.")
+        print(f"VC object stored under name: '{name}'.")
 
         # store the VC object dynamically
         store_VC = VirtualCircuit(
-            name=VC_name,
+            name=name,
             eq=eq,
             profiles=profiles,
             shape_matrix=shape_matrix,
             VCs_matrix=np.linalg.pinv(
                 shape_matrix
             ),  # "virtual circuits" are the pseudo-inverse of the shape matrix
-            targets=targets_new,
-            targets_val=self._targets_vec,
-            targets_options=targets_options,
-            non_standard_targets=non_standard_targets,
+            target_names=target_names,
             coils=coils,
+            target_calculator=target_calculator,
         )
-        setattr(self, VC_name, store_VC)
-
-        print("---")
-        print("Shape and virtual circuit matrices built.")
-        print(f"VC object stored under name: '{VC_name}'.")
+        setattr(self, name, store_VC)
 
     def apply_VC(
         self,
         eq,
         profiles,
         VC_object,
-        all_requested_target_shifts,
+        requested_target_shifts,
         verbose=False,
     ):
         """
@@ -676,11 +494,9 @@ def apply_VC(
             The profiles object upon which to apply the VCs.
         VC_object : an instance of the VirtualCircuit class
             Specifies the virtual circuit matrix and properties.
-        all_requested_targets_shift : list
+        requested_target_shifts : list
             List of floats containing the shifts in all of the relevant targets.
-            Same order as VC_object.targets.
-            Includes both standard and non-standard targets.
-            Functions to calculate non-standard targets are sourced from the VC_object.
+            Same order as VC_object.target_names.
         verbose: bool
             Display output (or not).
 
@@ -690,8 +506,6 @@ def apply_VC(
             Returns the equilibrium object after applying the shifted currents.
         object
             Returns the profiles object after applying the shifted currents.
-        list
-            List of strings containing all of targets of interest.
         np.array
             Array of new target values.
         np.array
@@ -699,16 +513,17 @@ def apply_VC(
         """
 
         # verify targets, coils, and shifts all match those used to generate VCs
-        assert len(all_requested_target_shifts) == len(
-            VC_object.targets_val
-        ), "The vector of requested shifts does not match the list of targets associated with the supplied VC_object!"
-        shifts = all_requested_target_shifts
+        if len(requested_target_shifts) != VC_object.VCs_matrix.shape[1]:
+            raise ValueError(
+                "The length of 'requested_target_shifts' does not match the list of targets "
+                "associated with the supplied VC object!"
+            )
 
         # calculate current shifts required using shape matrix (for stability)
         # uses least squares solver to solve S*dI = dT
         # where dT are the target shifts and dI the current shifts
         current_shifts = np.linalg.lstsq(
-            VC_object.shape_matrix, np.array(shifts), rcond=None
+            VC_object.shape_matrix, np.array(requested_target_shifts), rcond=None
         )[0]
 
         if verbose:
@@ -716,21 +531,20 @@ def apply_VC(
             print(f"{VC_object.coils} = {current_shifts}.")
 
         # re-solve static GS problem (to make sure it's solved already)
-        self.solver.forward_solve(
-            eq=eq,
-            profiles=profiles,
-            target_relative_tolerance=self.target_relative_tolerance,
-        )
+        try:
+            self.solver.forward_solve(
+                eq=eq,
+                profiles=profiles,
+                target_relative_tolerance=self.target_relative_tolerance,
+                suppress=True,
+            )
+        except AttributeError:
+            raise AttributeError("Solver not defined. Call define_solver() first.")
 
         # calculate the targets
-        _, old_target_values = self.calculate_targets(
-            eq,
-            VC_object.targets[
-                0 : len(VC_object.targets) - VC_object.len_non_standard_targets
-            ],
-            VC_object.targets_options,
-            VC_object.non_standard_targets,
-        )
+        # if not hasattr(self, "target_calculator"):
+        #     self.target_calculator = VC_object.target_calculator
+        old_target_values = VC_object.target_calculator(eq)
 
         # store copies of the eq and profile objects
         eq_new = eq.create_auxiliary_equilibrium()
@@ -744,24 +558,23 @@ def apply_VC(
         self.assign_currents(new_currents, VC_object.coils, eq=eq_new)
 
         # solve for the new equilibrium
-        self.solver.forward_solve(
-            eq_new,
-            profiles_new,
-            target_relative_tolerance=self.target_relative_tolerance,
-        )
+        try:
+            self.solver.forward_solve(
+                eq_new,
+                profiles_new,
+                target_relative_tolerance=self.target_relative_tolerance,
+                suppress=True,
+            )
+        except AttributeError:
+            raise AttributeError("Solver not defined. Call define_solver() first.")
 
         # calculate new target values and the difference vs. the old
-        target_names, new_target_values = self.calculate_targets(
-            eq_new,
-            VC_object.targets[
-                0 : len(VC_object.targets) - VC_object.len_non_standard_targets
-            ],
-            VC_object.targets_options,
-            VC_object.non_standard_targets,
-        )
+        new_target_values = VC_object.target_calculator(eq_new)
 
         if verbose:
             print(f"Targets shifts from VCs:")
-            print(f"{target_names} = {new_target_values - old_target_values}.")
+            print(
+                f"{VC_object.target_names} = {new_target_values - old_target_values}."
+            )
 
-        return eq_new, profiles_new, target_names, new_target_values, old_target_values
+        return eq_new, profiles_new, new_target_values, old_target_values