diff --git a/examples/3_Advanced/QA_reactorscale_fixed_orientations.py b/examples/3_Advanced/QA_reactorscale_fixed_orientations.py
index 23a52472b..efc5fb538 100644
--- a/examples/3_Advanced/QA_reactorscale_fixed_orientations.py
+++ b/examples/3_Advanced/QA_reactorscale_fixed_orientations.py
@@ -273,7 +273,7 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list):
 
 LENGTH_WEIGHT = Weight(0.001)
 LENGTH_TARGET = 130
-LINK_WEIGHT = 0.0
+LINK_WEIGHT = 1e3
 CC_THRESHOLD = 0.8
 CC_WEIGHT = 1e1
 CS_THRESHOLD = 1.5
diff --git a/examples/3_Advanced/QH_reactorscale_DA.py b/examples/3_Advanced/QH_reactorscale_DA.py
index c5ce880e7..77975721c 100644
--- a/examples/3_Advanced/QH_reactorscale_DA.py
+++ b/examples/3_Advanced/QH_reactorscale_DA.py
@@ -39,10 +39,10 @@
 
 # Initialize the boundary magnetic surface:
 range_param = "half period"
-nphi = 32
-ntheta = 32
-poff = 2.25
-coff = 3.5
+nphi = 64
+ntheta = 64
+poff = 2.0
+coff = 2.0
 s = SurfaceRZFourier.from_vmec_input(filename, range=range_param, nphi=nphi, ntheta=ntheta)
 s_inner = SurfaceRZFourier.from_vmec_input(filename, range=range_param, nphi=nphi * 4, ntheta=ntheta * 4)
 s_outer = SurfaceRZFourier.from_vmec_input(filename, range=range_param, nphi=nphi * 4, ntheta=ntheta * 4)
@@ -85,7 +85,7 @@ def initialize_coils_QH(TEST_DIR, s):
     from simsopt.geo import curves_to_vtk
 
     # generate planar TF coils
-    ncoils = 3
+    ncoils = 2
     R0 = s.get_rc(0, 0) * 1
     R1 = s.get_rc(1, 0) * 3
     order = 4
@@ -176,23 +176,50 @@ def initialize_coils_QH(TEST_DIR, s):
 
 ncoils = len(base_curves)
 print('Ncoils = ', ncoils)
+coil_normals = np.zeros((ncoils, 3))
+plasma_points = s.gamma().reshape(-1, 3)
+plasma_unitnormals = s.unitnormal().reshape(-1, 3)
+for i in range(ncoils):
+    point = (base_curves[i].get_dofs()[-3:])
+    dists = np.sum((point - plasma_points) ** 2, axis=-1)
+    min_ind = np.argmin(dists)
+    coil_normals[i, :] = plasma_unitnormals[min_ind, :]
+    # coil_normals[i, :] = (plasma_points[min_ind, :] - point)
+coil_normals = coil_normals / np.linalg.norm(coil_normals, axis=-1)[:, None]
+# alphas = np.arctan2(
+#                 -coil_normals[:, 1], 
+#                 np.sqrt(coil_normals[:, 0] ** 2 + coil_normals[:, 2] ** 2))
+# deltas = np.arcsin(coil_normals[:, 0] / \
+#                             np.sqrt(coil_normals[:, 0] ** 2 + coil_normals[:, 2] ** 2))
+alphas = np.arcsin(
+                -coil_normals[:, 1], 
+                )
+deltas = np.arctan2(coil_normals[:, 0], coil_normals[:, 2])
 for i in range(len(base_curves)):
-    # base_curves[i].set('x' + str(2 * order + 1), np.random.rand(1) - 0.5)
-    # base_curves[i].set('x' + str(2 * order + 2), np.random.rand(1) - 0.5)
-    # base_curves[i].set('x' + str(2 * order + 3), np.random.rand(1) - 0.5)
-    # base_curves[i].set('x' + str(2 * order + 4), np.random.rand(1) - 0.5)
+    alpha2 = alphas[i] / 2.0
+    delta2 = deltas[i] / 2.0
+    calpha2 = np.cos(alpha2)
+    salpha2 = np.sin(alpha2)
+    cdelta2 = np.cos(delta2)
+    sdelta2 = np.sin(delta2)
+    base_curves[i].set('x' + str(2 * order + 1), calpha2 * cdelta2)
+    base_curves[i].set('x' + str(2 * order + 2), salpha2 * cdelta2)
+    base_curves[i].set('x' + str(2 * order + 3), calpha2 * sdelta2)
+    base_curves[i].set('x' + str(2 * order + 4), -salpha2 * sdelta2)
+    # Fix orientations of each coil
+    base_curves[i].fix('x' + str(2 * order + 1))
+    base_curves[i].fix('x' + str(2 * order + 2))
+    base_curves[i].fix('x' + str(2 * order + 3))
+    base_curves[i].fix('x' + str(2 * order + 4))
 
     # Fix shape of each coil
     for j in range(2 * order + 1):
         base_curves[i].fix('x' + str(j))
     # Fix center points of each coil
-    # base_curves[i].fix('x' + str(2 * order + 5))
-    # base_curves[i].fix('x' + str(2 * order + 6))
-    # base_curves[i].fix('x' + str(2 * order + 7))
+    base_curves[i].fix('x' + str(2 * order + 5))
+    base_curves[i].fix('x' + str(2 * order + 6))
+    base_curves[i].fix('x' + str(2 * order + 7))
 base_currents = [Current(1e-1) * 2e7 for i in range(ncoils)]
-# Fix currents in each coil
-# for i in range(ncoils):
-#     base_currents[i].fix_all()
 
 coils = coils_via_symmetries(base_curves, base_currents, s.nfp, True)
 base_coils = coils[:ncoils]
@@ -237,25 +264,24 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list):
 
 LENGTH_WEIGHT = Weight(0.001)
 LENGTH_TARGET = 130
-LINK_WEIGHT = 0.0
+LINK_WEIGHT = 1e3
 CC_THRESHOLD = 0.8
 CC_WEIGHT = 1e1
 CS_THRESHOLD = 1.5
 CS_WEIGHT = 1e2
 # Weight for the Coil Coil forces term
-FORCE_WEIGHT = Weight(0.0) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons
+FORCE_WEIGHT = Weight(1e-22) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons
 FORCE_WEIGHT2 = Weight(0.0) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons
-TORQUE_WEIGHT = Weight(0.0) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons
-TORQUE_WEIGHT2 = Weight(0.0) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons
+TORQUE_WEIGHT = Weight(1e-24) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons
+TORQUE_WEIGHT2 = Weight(1e-24) # Forces are in Newtons, and typical values are ~10^5, 10^6 Newtons
 # Directory for output
-OUT_DIR = ("./QH_reactorscale/")
-    #        n{:d}_p{:.2e}_c{:.2e}_lw{:.2e}_lt{:.2e}_lkw{:.2e}" + \
-    # "_cct{:.2e}_ccw{:.2e}_cst{:.2e}_csw{:.2e}_fw{:.2e}_fww{:2e}_tw{:.2e}_tww{:2e}/").format(
-    #     ncoils, poff, coff, LENGTH_WEIGHT.value, LENGTH_TARGET, LINK_WEIGHT, 
-    #     CC_THRESHOLD, CC_WEIGHT, CS_THRESHOLD, CS_WEIGHT, FORCE_WEIGHT.value, 
-    #     FORCE_WEIGHT2.value,
-    #     TORQUE_WEIGHT.value,
-    #     TORQUE_WEIGHT2.value)
+OUT_DIR = ("./QH_reactorscale_n{:d}_p{:.2e}_c{:.2e}_lw{:.2e}_lt{:.2e}_lkw{:.2e}" + \
+    "_cct{:.2e}_ccw{:.2e}_cst{:.2e}_csw{:.2e}_fw{:.2e}_fww{:2e}_tw{:.2e}_tww{:2e}/").format(
+        ncoils, poff, coff, LENGTH_WEIGHT.value, LENGTH_TARGET, LINK_WEIGHT, 
+        CC_THRESHOLD, CC_WEIGHT, CS_THRESHOLD, CS_WEIGHT, FORCE_WEIGHT.value, 
+        FORCE_WEIGHT2.value,
+        TORQUE_WEIGHT.value,
+        TORQUE_WEIGHT2.value)
 if os.path.exists(OUT_DIR):
     shutil.rmtree(OUT_DIR)
 os.makedirs(OUT_DIR, exist_ok=True)
@@ -310,7 +336,7 @@ def pointData_forces_torques(coils, allcoils, aprimes, bprimes, nturns_list):
 # coil-coil and coil-plasma distances should be between all coils
 
 ### Jcc below removed the dipoles!
-Jccdist = CurveCurveDistance(curves_TF, CC_THRESHOLD, num_basecurves=len(coils_TF))
+Jccdist = CurveCurveDistance(curves + curves_TF, CC_THRESHOLD, num_basecurves=len(coils + coils_TF))
 Jcsdist = CurveSurfaceDistance(curves_TF, s, CS_THRESHOLD)
 
 # While the coil array is not moving around, they cannot
@@ -500,11 +526,11 @@ def fun(dofs):
 print('dJtorques time = ', t2 - t1, ' s')
 
 n_saves = 1
-MAXITER = 300
+MAXITER = 400
 for i in range(1, n_saves + 1):
     print('Iteration ' + str(i) + ' / ' + str(n_saves))
     res = minimize(fun, dofs, jac=True, method='L-BFGS-B', 
-        options={'maxiter': MAXITER, 'maxcor': 300}, tol=1e-15)
+        options={'maxiter': MAXITER, 'maxcor': 400}, tol=1e-15)
     # dofs = res.x
 
     dipole_currents = [c.current.get_value() for c in bs.coils]
@@ -550,6 +576,6 @@ def fun(dofs):
 
 t2 = time.time()
 print('Total time = ', t2 - t1)
-btot.save("biot_savart_optimized_QA.json")
+btot.save(OUT_DIR + "biot_savart_optimized_QH.json")
 print(OUT_DIR)