diff --git a/quartical/apps/plotter.py b/quartical/apps/plotter.py
index 35f18977..45577001 100644
--- a/quartical/apps/plotter.py
+++ b/quartical/apps/plotter.py
@@ -193,6 +193,9 @@ def to_plot_dict(xdsl, iter_attrs):
 
 
 def _plot(group, xds, args):
+    # get rid of question marks
+    qstrip = lambda x: x.replace('?', 'N/A')
+    group = tuple(map(qstrip, group))
 
     xds = xds.compute(scheduler="single-threaded")
 
diff --git a/quartical/calibration/constructor.py b/quartical/calibration/constructor.py
index ec6b6412..34462fac 100644
--- a/quartical/calibration/constructor.py
+++ b/quartical/calibration/constructor.py
@@ -5,7 +5,7 @@
 
 
 term_spec_tup = namedtuple("term_spec_tup", "name type shape pshape")
-aux_info_fields = ("SCAN_NUMBER", "FIELD_ID", "DATA_DESC_ID")
+log_info_fields = ("SCAN_NUMBER", "FIELD_ID", "DATA_DESC_ID")
 
 
 def construct_solver(
@@ -54,9 +54,9 @@ def construct_solver(
         corr_mode = data_xds.sizes["corr"]
 
         block_id_arr = get_block_id_arr(data_col)
-        aux_block_info = {
-            k: data_xds.attrs.get(k, "?") for k in aux_info_fields
-        }
+        data_xds_meta = data_xds.attrs.copy()
+        for k in log_info_fields:
+            data_xds_meta[k] = data_xds_meta.get(k, "?")
 
         # Grab the number of input chunks - doing this on the data should be
         # safe.
@@ -87,7 +87,7 @@ def construct_solver(
         )
         blocker.add_input("term_spec_list", spec_list, ("row", "chan"))
         blocker.add_input("corr_mode", corr_mode)
-        blocker.add_input("aux_block_info", aux_block_info)
+        blocker.add_input("data_xds_meta", data_xds_meta)
         blocker.add_input("solver_opts", solver_opts)
         blocker.add_input("chain", chain)
 
diff --git a/quartical/calibration/solver.py b/quartical/calibration/solver.py
index 620874c3..ddf3452f 100644
--- a/quartical/calibration/solver.py
+++ b/quartical/calibration/solver.py
@@ -54,7 +54,7 @@ def solver_wrapper(
     solver_opts,
     chain,
     block_id_arr,
-    aux_block_info,
+    data_xds_meta,
     corr_mode,
     **kwargs
 ):
@@ -108,11 +108,11 @@ def solver_wrapper(
         # Perform term specific setup e.g. init gains and params.
         if term.is_parameterized:
             gains, gain_flags, params, param_flags = term.init_term(
-                term_spec, ref_ant, ms_kwargs, term_kwargs
+                term_spec, ref_ant, ms_kwargs, term_kwargs, meta=data_xds_meta
             )
         else:
             gains, gain_flags = term.init_term(
-                term_spec, ref_ant, ms_kwargs, term_kwargs
+                term_spec, ref_ant, ms_kwargs, term_kwargs, meta=data_xds_meta
             )
             # Dummy arrays with standard dtypes - aids compilation.
             params = np.empty(term_pshape, dtype=np.float64)
@@ -190,6 +190,7 @@ def solver_wrapper(
     for ind, (term, iters) in enumerate(zip(cycle(chain), iter_recipe)):
 
         active_term = chain.index(term)
+        active_spec = term_spec_list[term_ind]
 
         ms_fields = term.ms_inputs._fields
         ms_inputs = term.ms_inputs(
@@ -219,13 +220,17 @@ def solver_wrapper(
             term.solve_per
         )
 
-        jhj, conv_iter, conv_perc = term.solver(
-            ms_inputs,
-            mapping_inputs,
-            chain_inputs,
-            meta_inputs,
-            corr_mode
-        )
+        if term.solver:
+            jhj, conv_iter, conv_perc = term.solver(
+                ms_inputs,
+                mapping_inputs,
+                chain_inputs,
+                meta_inputs,
+                corr_mode
+            )
+        else:
+            jhj = np.zeros(getattr(active_spec, "pshape", active_spec.shape))
+            conv_iter, conv_perc = 0, 1
 
         # If reweighting is enabled, do it when the epoch changes, except
         # for the final epoch - we don't reweight if we won't solve again.
@@ -269,7 +274,7 @@ def solver_wrapper(
         corr_mode
     )
 
-    log_chisq(presolve_chisq, postsolve_chisq, aux_block_info, block_id)
+    log_chisq(presolve_chisq, postsolve_chisq, data_xds_meta, block_id)
 
     results_dict["presolve_chisq"] = presolve_chisq
     results_dict["postsolve_chisq"] = postsolve_chisq
diff --git a/quartical/config/gain_schema.yaml b/quartical/config/gain_schema.yaml
index 3c3d0fc9..266c23b7 100644
--- a/quartical/config/gain_schema.yaml
+++ b/quartical/config/gain_schema.yaml
@@ -14,7 +14,9 @@ gain:
       - rotation_measure
       - rotation
       - crosshand_phase
+      - crosshand_phase_null_v
       - leakage
+      - parallactic_angle
     info:
       Type of gain to solve for.
 
diff --git a/quartical/data_handling/angles.py b/quartical/data_handling/angles.py
index 62ca7c34..aa05b1be 100644
--- a/quartical/data_handling/angles.py
+++ b/quartical/data_handling/angles.py
@@ -19,6 +19,46 @@
 _thread_local = threading.local()
 
 
+def assign_parangle_data(ms_path, data_xds_list):
+
+    anttab = xds_from_storage_table(ms_path + "::ANTENNA")[0]
+    feedtab = xds_from_storage_table(ms_path + "::FEED")[0]
+    fieldtab = xds_from_storage_table(ms_path + "::FIELD")[0]
+
+    # We do the following eagerly to reduce graph complexity.
+    feeds = feedtab.POLARIZATION_TYPE.values
+    unique_feeds = np.unique(feeds)
+
+    if np.all([feed in "XxYy" for feed in unique_feeds]):
+        feed_type = "linear"
+    elif np.all([feed in "LlRr" for feed in unique_feeds]):
+        feed_type = "circular"
+    else:
+        raise ValueError("Unsupported feed type/configuration.")
+
+    phase_dirs = fieldtab.PHASE_DIR.values
+
+    updated_data_xds_list = []
+    for xds in data_xds_list:
+        xds = xds.assign(
+            {
+                "RECEPTOR_ANGLE": (
+                    ("ant", "feed"), clone(feedtab.RECEPTOR_ANGLE.data)
+                ),
+                "POSITION": (
+                    ("ant", "xyz"),
+                    clone(anttab.POSITION.data)
+                )
+            }
+        )
+        xds.attrs["FEED_TYPE"] = feed_type
+        xds.attrs["FIELD_CENTRE"] = tuple(phase_dirs[xds.FIELD_ID, 0])
+
+        updated_data_xds_list.append(xds)
+
+    return updated_data_xds_list
+
+
 def make_parangle_xds_list(ms_path, data_xds_list):
     """Create a list of xarray.Datasets containing the parallactic angles."""
 
@@ -266,7 +306,7 @@ def nb_apply_parangle_rot(data_col, parangles, utime_ind, ant1_col, ant2_col,
     v1_imul_v2 = factories.v1_imul_v2_factory(corr_mode)
     v1_imul_v2ct = factories.v1_imul_v2ct_factory(corr_mode)
     valloc = factories.valloc_factory(corr_mode)
-    rotmat = rotation_factory(corr_mode, feed_type)
+    rotmat = factories.rotation_factory(corr_mode, feed_type)
 
     def impl(data_col, parangles, utime_ind, ant1_col, ant2_col,
              corr_mode, feed_type):
@@ -299,43 +339,3 @@ def impl(data_col, parangles, utime_ind, ant1_col, ant2_col,
         return data_col
 
     return impl
-
-
-def rotation_factory(corr_mode, feed_type):
-
-    if feed_type.literal_value == "circular":
-        if corr_mode.literal_value == 4:
-            def impl(rot0, rot1, out):
-                out[0] = np.exp(-1j*rot0)
-                out[1] = 0
-                out[2] = 0
-                out[3] = np.exp(1j*rot1)
-        elif corr_mode.literal_value == 2:  # TODO: Is this sensible?
-            def impl(rot0, rot1, out):
-                out[0] = np.exp(-1j*rot0)
-                out[1] = np.exp(1j*rot1)
-        elif corr_mode.literal_value == 1:  # TODO: Is this sensible?
-            def impl(rot0, rot1, out):
-                out[0] = np.exp(-1j*rot0)
-        else:
-            raise ValueError("Unsupported number of correlations.")
-    elif feed_type.literal_value == "linear":
-        if corr_mode.literal_value == 4:
-            def impl(rot0, rot1, out):
-                out[0] = np.cos(rot0)
-                out[1] = np.sin(rot0)
-                out[2] = -np.sin(rot1)
-                out[3] = np.cos(rot1)
-        elif corr_mode.literal_value == 2:  # TODO: Is this sensible?
-            def impl(rot0, rot1, out):
-                out[0] = np.cos(rot0)
-                out[1] = np.cos(rot1)
-        elif corr_mode.literal_value == 1:  # TODO: Is this sensible?
-            def impl(rot0, rot1, out):
-                out[0] = np.cos(rot0)
-        else:
-            raise ValueError("Unsupported number of correlations.")
-    else:
-        raise ValueError("Unsupported feed type.")
-
-    return factories.qcjit(impl)
diff --git a/quartical/data_handling/ms_handler.py b/quartical/data_handling/ms_handler.py
index 3f5b38cf..13920fcb 100644
--- a/quartical/data_handling/ms_handler.py
+++ b/quartical/data_handling/ms_handler.py
@@ -15,6 +15,13 @@
 from quartical.data_handling.selection import filter_xds_list
 from quartical.data_handling.angles import apply_parangles
 
+DASKMS_ATTRS = {
+    "__daskms_partition_schema__",
+    "SCAN_NUMBER",
+    "FIELD_ID",
+    "DATA_DESC_ID"
+}
+
 
 def read_xds_list(model_columns, ms_opts):
     """Reads a measurement set and generates a list of xarray data sets.
@@ -293,14 +300,18 @@ def write_xds_list(xds_list, ref_xds_list, ms_path, output_opts):
 
     logger.info("Outputs will be written to {}.", ", ".join(output_cols))
 
+    # Select only the output columns to simplify datasets.
+    xds_list = [xds[list(output_cols)] for xds in xds_list]
+
+    # Remove all coords bar ROWID so that they do not get written.
+    xds_list = [
+        xds.drop_vars(set(xds.coords.keys()) - {"ROWID"}, errors='ignore')
+        for xds in xds_list
+    ]
+
     # Remove attrs added by QuartiCal so that they do not get written.
     for xds in xds_list:
-        xds.attrs.pop("UTIME_CHUNKS", None)
-        xds.attrs.pop("FIELD_NAME", None)
-
-    # Remove coords added by QuartiCal so that they do not get written.
-    xds_list = [xds.drop_vars(["chan", "corr"], errors='ignore')
-                for xds in xds_list]
+        xds.attrs = {k: v for k, v in xds.attrs.items() if k in DASKMS_ATTRS}
 
     with warnings.catch_warnings():  # We anticipate spurious warnings.
         warnings.simplefilter("ignore")
diff --git a/quartical/executor.py b/quartical/executor.py
index b79a9fe8..704ecdaf 100644
--- a/quartical/executor.py
+++ b/quartical/executor.py
@@ -13,7 +13,8 @@
                                                 preprocess_xds_list,
                                                 postprocess_xds_list)
 from quartical.data_handling.model_handler import add_model_graph
-from quartical.data_handling.angles import make_parangle_xds_list
+from quartical.data_handling.angles import (make_parangle_xds_list,
+                                            assign_parangle_data)
 from quartical.calibration.calibrate import add_calibration_graph
 from quartical.statistics.statistics import make_stats_xds_list
 from quartical.statistics.logging import log_summary_stats
@@ -110,6 +111,7 @@ def _execute(exitstack):
 
     # Preprocess the xds_list - initialise some values and fix bad data.
     data_xds_list = preprocess_xds_list(data_xds_list, ms_opts)
+    data_xds_list = assign_parangle_data(ms_opts.path, data_xds_list)
 
     # Make a list of datasets containing the parallactic angles as these
     # can be expensive to compute and may be used several times. NOTE: At
diff --git a/quartical/gains/__init__.py b/quartical/gains/__init__.py
index e799443c..e652419d 100644
--- a/quartical/gains/__init__.py
+++ b/quartical/gains/__init__.py
@@ -6,9 +6,10 @@
 from quartical.gains.tec_and_offset import TecAndOffset
 from quartical.gains.rotation import Rotation
 from quartical.gains.rotation_measure import RotationMeasure
-from quartical.gains.crosshand_phase import CrosshandPhase
+from quartical.gains.crosshand_phase import CrosshandPhase, CrosshandPhaseNullV
 from quartical.gains.leakage import Leakage
 from quartical.gains.delay_and_tec import DelayAndTec
+from quartical.gains.parallactic_angle import ParallacticAngle
 
 
 TERM_TYPES = {
@@ -22,6 +23,8 @@
     "rotation": Rotation,
     "rotation_measure": RotationMeasure,
     "crosshand_phase": CrosshandPhase,
+    "crosshand_phase_null_v": CrosshandPhaseNullV,
     "leakage": Leakage,
-    "delay_and_tec": DelayAndTec
+    "delay_and_tec": DelayAndTec,
+    "parallactic_angle": ParallacticAngle
 }
diff --git a/quartical/gains/amplitude/__init__.py b/quartical/gains/amplitude/__init__.py
index 5e04ba09..a2d2285a 100644
--- a/quartical/gains/amplitude/__init__.py
+++ b/quartical/gains/amplitude/__init__.py
@@ -38,7 +38,7 @@ def make_param_names(cls, correlations):
 
         return [f"amplitude_{c}" for c in param_corr]
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
diff --git a/quartical/gains/crosshand_phase/__init__.py b/quartical/gains/crosshand_phase/__init__.py
index 80fd8e9c..f9296355 100644
--- a/quartical/gains/crosshand_phase/__init__.py
+++ b/quartical/gains/crosshand_phase/__init__.py
@@ -5,6 +5,9 @@
     crosshand_phase_solver,
     crosshand_params_to_gains
 )
+from quartical.gains.crosshand_phase.null_v_kernel import (
+    null_v_crosshand_phase_solver
+)
 from quartical.gains.general.flagging import (
     apply_gain_flags_to_gains,
     apply_param_flags_to_params
@@ -39,7 +42,7 @@ def make_param_names(cls, correlations):
 
         return [f"crosshand_phase_{c}" for c in param_corr]
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
@@ -54,3 +57,8 @@ def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
         apply_gain_flags_to_gains(gain_flags, gains)
 
         return gains, gain_flags, params, param_flags
+
+
+class CrosshandPhaseNullV(CrosshandPhase):
+
+    solver = staticmethod(null_v_crosshand_phase_solver)
\ No newline at end of file
diff --git a/quartical/gains/crosshand_phase/null_v_kernel.py b/quartical/gains/crosshand_phase/null_v_kernel.py
new file mode 100644
index 00000000..33fb4f78
--- /dev/null
+++ b/quartical/gains/crosshand_phase/null_v_kernel.py
@@ -0,0 +1,643 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+from numba import prange, njit
+from numba.typed import List
+from numba.extending import overload
+from quartical.utils.numba import (coerce_literal,
+                                   JIT_OPTIONS,
+                                   PARALLEL_JIT_OPTIONS)
+from quartical.gains.general.generics import (native_intermediaries,
+                                              upsampled_itermediaries,
+                                              per_array_jhj_jhr,
+                                              resample_solints,
+                                              downsample_jhj_jhr,
+                                              invert_gains)
+from quartical.gains.general.flagging import (flag_intermediaries,
+                                              update_gain_flags,
+                                              finalize_gain_flags,
+                                              apply_gain_flags_to_flag_col,
+                                              update_param_flags)
+from quartical.gains.general.convenience import (get_row,
+                                                 get_extents)
+import quartical.gains.general.factories as factories
+from quartical.gains.general.inversion import (invert_factory,
+                                               inversion_buffer_factory)
+
+
+def get_identity_params(corr_mode):
+
+    if corr_mode.literal_value == 4:
+        return np.zeros((1,), dtype=np.float64)
+    else:
+        raise ValueError("Unsupported number of correlations.")
+
+
+@njit(**JIT_OPTIONS)
+def null_v_crosshand_phase_solver(
+    ms_inputs,
+    mapping_inputs,
+    chain_inputs,
+    meta_inputs,
+    corr_mode
+):
+    return null_v_crosshand_phase_solver_impl(
+        ms_inputs,
+        mapping_inputs,
+        chain_inputs,
+        meta_inputs,
+        corr_mode
+    )
+
+
+def null_v_crosshand_phase_solver_impl(
+    ms_inputs,
+    mapping_inputs,
+    chain_inputs,
+    meta_inputs,
+    corr_mode
+):
+    raise NotImplementedError
+
+
+@overload(null_v_crosshand_phase_solver_impl, jit_options=JIT_OPTIONS)
+def nb_null_v_crosshand_phase_solver_impl(
+    ms_inputs,
+    mapping_inputs,
+    chain_inputs,
+    meta_inputs,
+    corr_mode
+):
+
+    coerce_literal(nb_null_v_crosshand_phase_solver_impl, ["corr_mode"])
+
+    identity_params = get_identity_params(corr_mode)
+
+    def impl(
+        ms_inputs,
+        mapping_inputs,
+        chain_inputs,
+        meta_inputs,
+        corr_mode
+    ):
+
+        gains = chain_inputs.gains
+        gain_flags = chain_inputs.gain_flags
+
+        inverse_gains = List()
+        for gain_term in gains:
+            inverse_gains.append(np.empty_like(gain_term))
+        invert_gains(gains, inverse_gains, corr_mode)
+
+        active_term = meta_inputs.active_term
+        max_iter = meta_inputs.iters
+        solve_per = meta_inputs.solve_per
+        dd_term = meta_inputs.dd_term
+        n_thread = meta_inputs.threads
+
+        active_gain = gains[active_term]
+        active_gain_flags = gain_flags[active_term]
+        active_params = chain_inputs.params[active_term]
+
+        # Set up some intemediaries used for flagging.
+        km1_gain = active_gain.copy()
+        km1_abs2_diffs = np.zeros_like(active_gain_flags, dtype=np.float64)
+        abs2_diffs_trend = np.zeros_like(active_gain_flags, dtype=np.float64)
+        flag_imdry = \
+            flag_intermediaries(km1_gain, km1_abs2_diffs, abs2_diffs_trend)
+
+        # Set up some intemediaries used for solving.
+        real_dtype = active_gain.real.dtype
+        param_shape = active_params.shape
+
+        active_t_map_g = mapping_inputs.time_maps[active_term]
+        active_f_map_g = mapping_inputs.freq_maps[active_term]
+
+        # Create more work to do in paralllel when needed, else no-op.
+        resampler = resample_solints(active_t_map_g, param_shape, n_thread)
+
+        # Determine the starts and stops of the rows and channels associated
+        # with each solution interval.
+        extents = get_extents(resampler.upsample_t_map, active_f_map_g)
+
+        upsample_shape = resampler.upsample_shape
+        upsampled_jhj = np.empty(upsample_shape + (upsample_shape[-1],),
+                                 dtype=real_dtype)
+        upsampled_jhr = np.empty(upsample_shape, dtype=real_dtype)
+        jhj = upsampled_jhj[:param_shape[0]]
+        jhr = upsampled_jhr[:param_shape[0]]
+        update = np.zeros(param_shape, dtype=real_dtype)
+
+        upsampled_imdry = upsampled_itermediaries(upsampled_jhj, upsampled_jhr)
+        native_imdry = native_intermediaries(jhj, jhr, update)
+
+        for loop_idx in range(max_iter or 1):
+
+            compute_jhj_jhr(
+                inverse_gains,
+                ms_inputs,
+                mapping_inputs,
+                chain_inputs,
+                meta_inputs,
+                upsampled_imdry,
+                extents,
+                corr_mode
+            )
+
+            if resampler.active:
+                downsample_jhj_jhr(upsampled_imdry, resampler.downsample_t_map)
+
+            if solve_per == "array":
+                per_array_jhj_jhr(native_imdry)
+
+            if not max_iter:  # Non-solvable term, we just want jhj.
+                conv_perc = 0  # Didn't converge.
+                loop_idx = -1  # Did zero iterations.
+                break
+
+            compute_update(native_imdry, corr_mode)
+
+            finalize_update(
+                chain_inputs,
+                meta_inputs,
+                native_imdry,
+                loop_idx,
+                corr_mode
+            )
+
+            # The parameters/gains are correct but we are solving for the
+            # inverse so we update the inverse term here.
+            inverse_gains[active_term][:] = active_gain.conj()
+
+            # Check for gain convergence. Produced as a side effect of
+            # flagging. The converged percentage is based on unflagged
+            # intervals.
+            conv_perc = update_gain_flags(
+                chain_inputs,
+                meta_inputs,
+                flag_imdry,
+                loop_idx,
+                corr_mode,
+                numbness=1e9
+            )
+
+            # Propagate gain flags to parameter flags.
+            update_param_flags(
+                mapping_inputs,
+                chain_inputs,
+                meta_inputs,
+                identity_params
+            )
+
+            if conv_perc >= meta_inputs.stop_frac:
+                break
+
+        # NOTE: Removes soft flags and flags points which have bad trends.
+        finalize_gain_flags(
+            chain_inputs,
+            meta_inputs,
+            flag_imdry,
+            corr_mode
+        )
+
+        # Call this one last time to ensure points flagged by finialize are
+        # propagated (in the DI case).
+        if not dd_term:
+            apply_gain_flags_to_flag_col(
+                ms_inputs,
+                mapping_inputs,
+                chain_inputs,
+                meta_inputs
+            )
+
+        return native_imdry.jhj, loop_idx + 1, conv_perc
+
+    return impl
+
+
+def compute_jhj_jhr(
+    inverse_gains,
+    ms_inputs,
+    mapping_inputs,
+    chain_inputs,
+    meta_inputs,
+    upsampled_imdry,
+    extents,
+    corr_mode
+):
+    return NotImplementedError
+
+
+@overload(compute_jhj_jhr, jit_options=PARALLEL_JIT_OPTIONS)
+def nb_compute_jhj_jhr(
+    inverse_gains,
+    ms_inputs,
+    mapping_inputs,
+    chain_inputs,
+    meta_inputs,
+    upsampled_imdry,
+    extents,
+    corr_mode
+):
+
+    coerce_literal(nb_compute_jhj_jhr, ["corr_mode"])
+
+    # We want to dispatch based on this field so we need its type.
+    row_weights_idx = ms_inputs.fields.index('ROW_WEIGHTS')
+    row_weights_type = ms_inputs[row_weights_idx]
+
+    imul_rweight = factories.imul_rweight_factory(corr_mode, row_weights_type)
+    v1_imul_v2 = factories.v1_imul_v2_factory(corr_mode)
+    v1_imul_v2ct = factories.v1_imul_v2ct_factory(corr_mode)
+    v1ct_imul_v2 = factories.v1ct_imul_v2_factory(corr_mode)
+    absv1_idiv_absv2 = factories.absv1_idiv_absv2_factory(corr_mode)
+    iunpack = factories.iunpack_factory(corr_mode)
+    iunpackct = factories.iunpackct_factory(corr_mode)
+    imul = factories.imul_factory(corr_mode)
+    iadd = factories.iadd_factory(corr_mode)
+    isub = factories.isub_factory(corr_mode)
+    valloc = factories.valloc_factory(corr_mode)
+    make_loop_vars = factories.loop_var_factory(corr_mode)
+    set_identity = factories.set_identity_factory(corr_mode)
+    compute_jhwj_jhwr_elem = compute_jhwj_jhwr_elem_factory(corr_mode)
+
+    def impl(
+        inverse_gains,
+        ms_inputs,
+        mapping_inputs,
+        chain_inputs,
+        meta_inputs,
+        upsampled_imdry,
+        extents,
+        corr_mode
+    ):
+
+        data = ms_inputs.DATA
+        model = ms_inputs.MODEL_DATA
+        weights = ms_inputs.WEIGHT
+        flags = ms_inputs.FLAG
+        antenna1 = ms_inputs.ANTENNA1
+        antenna2 = ms_inputs.ANTENNA2
+        row_map = ms_inputs.ROW_MAP
+        row_weights = ms_inputs.ROW_WEIGHTS
+
+        # Reverse all the (invernse) gains and mappings.
+        time_maps = mapping_inputs.time_maps[::-1]
+        freq_maps = mapping_inputs.freq_maps[::-1]
+        dir_maps = mapping_inputs.dir_maps[::-1]
+
+        gains = inverse_gains[::-1]
+
+        # Active term in the REVERSED chain.
+        active_term = len(gains) - meta_inputs.active_term - 1
+
+        jhj = upsampled_imdry.jhj
+        jhr = upsampled_imdry.jhr
+
+        _, n_chan, n_dir, n_corr = model.shape
+
+        jhj[:] = 0
+        jhr[:] = 0
+
+        n_tint, n_fint, n_ant, n_gdir, n_param = jhr.shape
+        n_int = n_tint*n_fint
+
+        complex_dtype = gains[active_term].dtype
+        weight_dtype = weights.dtype
+
+        n_gains = len(gains)
+
+        row_starts = extents.row_starts
+        row_stops = extents.row_stops
+        chan_starts = extents.chan_starts
+        chan_stops = extents.chan_stops
+
+        # Determine loop variables based on where we are in the chain.
+        # gt means greater than (n>j) and lt means less than (n<j).
+        all_terms, gt_active, lt_active = make_loop_vars(n_gains, active_term)
+
+        # Parallel over all solution intervals.
+        for i in prange(n_int):
+
+            ti = i//n_fint
+            fi = i - ti*n_fint
+
+            rs = row_starts[ti]
+            re = row_stops[ti]
+            fs = chan_starts[fi]
+            fe = chan_stops[fi]
+
+            rop_pq = valloc(complex_dtype)  # Right-multiply operator for pq.
+            rop_qp = valloc(complex_dtype)  # Right-multiply operator for qp.
+            lop_pq = valloc(complex_dtype)  # Left-multiply operator for pq.
+            lop_qp = valloc(complex_dtype)  # Left-multiply operator for qp.
+
+            w = valloc(weight_dtype)
+            r_pq = valloc(complex_dtype)
+            wr_pq = valloc(complex_dtype)
+            wr_qp = valloc(complex_dtype)
+            v_pqd = valloc(complex_dtype)
+            v_pq = valloc(complex_dtype)
+
+            gains_p = valloc(complex_dtype, leading_dims=(n_gains,))
+            gains_q = valloc(complex_dtype, leading_dims=(n_gains,))
+
+            lop_pq_arr = valloc(complex_dtype, leading_dims=(n_gdir,))
+            rop_pq_arr = valloc(complex_dtype, leading_dims=(n_gdir,))
+            lop_qp_arr = valloc(complex_dtype, leading_dims=(n_gdir,))
+            rop_qp_arr = valloc(complex_dtype, leading_dims=(n_gdir,))
+
+            norm_factors = valloc(complex_dtype)
+
+            jhr_tifi = jhr[ti, fi]
+            jhj_tifi = jhj[ti, fi]
+
+            for row_ind in range(rs, re):
+
+                row = get_row(row_ind, row_map)
+                a1_m, a2_m = antenna1[row], antenna2[row]
+
+                for f in range(fs, fe):
+
+                    if flags[row, f]:  # Skip flagged data points.
+                        continue
+
+                    # Apply row weights in the BDA case, otherwise a no-op.
+                    imul_rweight(weights[row, f], w, row_weights, row_ind)
+
+                    lop_pq_arr[:] = 0
+                    rop_pq_arr[:] = 0
+                    lop_qp_arr[:] = 0
+                    rop_qp_arr[:] = 0
+                    v_pq[:] = 0
+                    r_pq[:] = 0
+
+                    for d in range(n_dir):
+
+                        set_identity(lop_pq)
+                        set_identity(lop_qp)
+
+                        # Construct a small contiguous gain array. This makes
+                        # the single term case fractionally slower.
+                        for gi in range(n_gains):
+                            d_m = dir_maps[gi][d]  # Broadcast dir.
+                            t_m = time_maps[gi][row_ind]
+                            f_m = freq_maps[gi][f]
+
+                            gain = gains[gi][t_m, f_m]
+
+                            iunpack(gains_p[gi], gain[a1_m, d_m])
+                            iunpack(gains_q[gi], gain[a2_m, d_m])
+
+                        m = data[row, f]
+                        iunpack(rop_qp, m)
+                        iunpackct(rop_pq, rop_qp)
+
+                        for g in all_terms:
+
+                            g_q = gains_q[g]
+                            v1_imul_v2(g_q, rop_pq, rop_pq)
+
+                            g_p = gains_p[g]
+                            v1_imul_v2(g_p, rop_qp, rop_qp)
+
+                        for g in gt_active:
+
+                            g_p = gains_p[g]
+                            v1_imul_v2ct(rop_pq, g_p, rop_pq)
+
+                            g_q = gains_q[g]
+                            v1_imul_v2ct(rop_qp, g_q, rop_qp)
+
+                        for g in lt_active:
+
+                            g_p = gains_p[g]
+                            v1ct_imul_v2(g_p, lop_pq, lop_pq)
+
+                            g_q = gains_q[g]
+                            v1ct_imul_v2(g_q, lop_qp, lop_qp)
+
+                        out_d = dir_maps[active_term][d]
+
+                        iunpack(lop_pq_arr[out_d], lop_pq)
+                        iadd(rop_pq_arr[out_d], rop_pq)
+
+                        iunpack(lop_qp_arr[out_d], lop_qp)
+                        iadd(rop_qp_arr[out_d], rop_qp)
+
+                        v1ct_imul_v2(lop_pq, gains_p[active_term], v_pqd)
+                        v1_imul_v2ct(v_pqd, rop_pq, v_pqd)
+                        iadd(v_pq, v_pqd)
+
+                    isub(r_pq, v_pq)
+
+                    for d in range(n_gdir):
+
+                        iunpack(wr_pq, r_pq)
+                        iunpackct(wr_qp, wr_pq)
+
+                        lop_pq_d = lop_pq_arr[d]
+                        rop_pq_d = rop_pq_arr[d]
+
+                        compute_jhwj_jhwr_elem(lop_pq_d,
+                                               rop_pq_d,
+                                               w,
+                                               norm_factors,
+                                               gains_p[active_term],
+                                               wr_pq,
+                                               jhr_tifi[a1_m, d],
+                                               jhj_tifi[a1_m, d])
+
+                        lop_qp_d = lop_qp_arr[d]
+                        rop_qp_d = rop_qp_arr[d]
+
+                        compute_jhwj_jhwr_elem(lop_qp_d,
+                                               rop_qp_d,
+                                               w,
+                                               norm_factors,
+                                               gains_q[active_term],
+                                               wr_qp,
+                                               jhr_tifi[a2_m, d],
+                                               jhj_tifi[a2_m, d])
+        return
+    return impl
+
+
+def compute_update(native_imdry, corr_mode):
+    raise NotImplementedError
+
+
+@overload(compute_update, jit_options=PARALLEL_JIT_OPTIONS)
+def nb_compute_update(native_imdry, corr_mode):
+
+    coerce_literal(nb_compute_update, ["corr_mode"])
+
+    # We want to dispatch based on this field so we need its type.
+    jhj = native_imdry[native_imdry.fields.index('jhj')]
+
+    generalised = jhj.ndim == 6
+    inversion_buffer = inversion_buffer_factory(generalised=generalised)
+    invert = invert_factory(corr_mode, generalised=generalised)
+
+    def impl(native_imdry, corr_mode):
+
+        jhj = native_imdry.jhj
+        jhr = native_imdry.jhr
+        update = native_imdry.update
+
+        n_tint, n_fint, n_ant, n_dir, n_param = jhr.shape
+
+        n_int = n_tint * n_fint
+
+        result_dtype = jhr.dtype
+
+        for i in prange(n_int):
+
+            t = i // n_fint
+            f = i - t * n_fint
+
+            buffers = inversion_buffer(n_param, result_dtype)
+
+            for a in range(n_ant):
+                for d in range(n_dir):
+
+                    invert(jhj[t, f, a, d],
+                           jhr[t, f, a, d],
+                           update[t, f, a, d],
+                           buffers)
+
+    return impl
+
+
+def finalize_update(
+    chain_inputs,
+    meta_inputs,
+    native_imdry,
+    loop_idx,
+    corr_mode
+):
+    raise NotImplementedError
+
+
+@overload(finalize_update, jit_options=JIT_OPTIONS)
+def nb_finalize_update(
+    chain_inputs,
+    meta_inputs,
+    native_imdry,
+    loop_idx,
+    corr_mode
+):
+
+    coerce_literal(nb_finalize_update, ["corr_mode"])
+
+    set_identity = factories.set_identity_factory(corr_mode)
+    param_to_gain = param_to_gain_factory(corr_mode)
+
+    def impl(
+        chain_inputs,
+        meta_inputs,
+        native_imdry,
+        loop_idx,
+        corr_mode
+    ):
+
+        dd_term = meta_inputs.dd_term
+        active_term = meta_inputs.active_term
+        pinned_directions = meta_inputs.pinned_directions
+
+        gain = chain_inputs.gains[active_term]
+        gain_flags = chain_inputs.gain_flags[active_term]
+        params = chain_inputs.params[active_term]
+
+        update = native_imdry.update
+
+        n_tint, n_fint, n_ant, n_dir, n_corr = gain.shape
+
+        if dd_term:
+            dir_loop = [d for d in range(n_dir) if d not in pinned_directions]
+        else:
+            dir_loop = [d for d in range(n_dir)]
+
+        for ti in range(n_tint):
+            for fi in range(n_fint):
+                for a in range(n_ant):
+                    for d in dir_loop:
+
+                        p = params[ti, fi, a, d]
+                        g = gain[ti, fi, a, d]
+                        fl = gain_flags[ti, fi, a, d]
+                        upd = update[ti, fi, a, d]
+
+                        if fl == 1:
+                            p[:] = 0
+                            set_identity(g)
+                        else:
+                            # NOTE: Halving the update absolutely required.
+                            p -= 0.5*upd  # Flip sign for non-inverse solution.
+                            param_to_gain(p, g)
+
+    return impl
+
+
+def param_to_gain_factory(corr_mode):
+
+    if corr_mode.literal_value == 4:
+        def impl(params, gain):
+            gain[0] = np.exp(1j*params[0])
+    else:
+        raise ValueError("Crosshand phase can only be solved for with four "
+                         "correlation data.")
+
+    return factories.qcjit(impl)
+
+
+def compute_jhwj_jhwr_elem_factory(corr_mode):
+
+    v1_imul_v2 = factories.v1_imul_v2_factory(corr_mode)
+    unpack = factories.unpack_factory(corr_mode)
+    unpackc = factories.unpackc_factory(corr_mode)
+    iabsdivsq = factories.iabsdivsq_factory(corr_mode)
+    imul = factories.imul_factory(corr_mode)
+
+    if corr_mode.literal_value == 4:
+        def impl(lop, rop, w, normf, gain, res, jhr, jhj):
+
+            _, r_1, r_2, _ = unpack(res)  # NOTE: XX, XY, YX, YY
+            v_res = -0.5j * r_1 + 0.5j * r_2
+
+            res[0] = 0
+            res[1] = 0.5j * v_res
+            res[2] = -0.5j * v_res
+            res[3] = 0
+
+            # Accumulate an element of jhwr.
+            v1_imul_v2(res, rop, res)
+            v1_imul_v2(lop, res, res)
+
+            r_0, _, _, _ = unpack(res)  # NOTE: XX, XY, YX, YY
+
+            gc_0, _, _, _ = unpackc(gain)
+
+            drv_00 = -1j*gc_0
+
+            upd_00 = (drv_00*r_0).real
+
+            jhr[0] += upd_00
+
+            lop_00, lop_01, _, _ = unpack(lop)
+            rop_00, _, rop_01, _ = unpack(rop)  # "Transpose"
+
+            jh_01 = lop_00 * rop_01
+            jh_02 = lop_01 * rop_00
+
+            jh_v = 0.5j * jh_01 - 0.5j * jh_02
+            j_v = jh_v.conjugate()
+
+            jhj_v = jh_v * j_v
+
+            jhj[0, 0] += jhj_v.real
+
+    else:
+        raise ValueError("Crosshand phase can only be solved for with four "
+                         "correlation data.")
+
+    return factories.qcjit(impl)
diff --git a/quartical/gains/delay/__init__.py b/quartical/gains/delay/__init__.py
index 289575d3..2e4a4ff5 100644
--- a/quartical/gains/delay/__init__.py
+++ b/quartical/gains/delay/__init__.py
@@ -52,7 +52,7 @@ def make_param_names(cls, correlations):
 
         return [n.format(c) for c in param_corr for n in template]
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
diff --git a/quartical/gains/delay_and_offset/__init__.py b/quartical/gains/delay_and_offset/__init__.py
index de88830a..0b74d3a2 100644
--- a/quartical/gains/delay_and_offset/__init__.py
+++ b/quartical/gains/delay_and_offset/__init__.py
@@ -55,7 +55,7 @@ def make_param_names(cls, correlations):
 
         return [n.format(c) for c in param_corr for n in template]
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
diff --git a/quartical/gains/delay_and_tec/__init__.py b/quartical/gains/delay_and_tec/__init__.py
index 9e7d5720..7d58e80a 100644
--- a/quartical/gains/delay_and_tec/__init__.py
+++ b/quartical/gains/delay_and_tec/__init__.py
@@ -54,7 +54,7 @@ def make_param_names(cls, correlations):
 
         return [n.format(c) for c in param_corr for n in template]
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
diff --git a/quartical/gains/gain.py b/quartical/gains/gain.py
index ca1165fd..0b18d4f2 100644
--- a/quartical/gains/gain.py
+++ b/quartical/gains/gain.py
@@ -344,7 +344,7 @@ def _make_dir_map(cls, n_dir, direction_dependent):
 
         return dir_map
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         (_, _, gain_shape, _) = term_spec
diff --git a/quartical/gains/general/factories.py b/quartical/gains/general/factories.py
index 469854a0..a73bd205 100644
--- a/quartical/gains/general/factories.py
+++ b/quartical/gains/general/factories.py
@@ -798,3 +798,43 @@ def impl(a, b, out):
         out[3, 3] = a11 * b11
 
     return qcjit(impl)
+
+
+def rotation_factory(corr_mode, feed_type):
+
+    if feed_type.literal_value == "circular":
+        if corr_mode.literal_value == 4:
+            def impl(rot0, rot1, out):
+                out[0] = np.exp(-1j*rot0)
+                out[1] = 0
+                out[2] = 0
+                out[3] = np.exp(1j*rot1)
+        elif corr_mode.literal_value == 2:  # TODO: Is this sensible?
+            def impl(rot0, rot1, out):
+                out[0] = np.exp(-1j*rot0)
+                out[1] = np.exp(1j*rot1)
+        elif corr_mode.literal_value == 1:  # TODO: Is this sensible?
+            def impl(rot0, rot1, out):
+                out[0] = np.exp(-1j*rot0)
+        else:
+            raise ValueError("Unsupported number of correlations.")
+    elif feed_type.literal_value == "linear":
+        if corr_mode.literal_value == 4:
+            def impl(rot0, rot1, out):
+                out[0] = np.cos(rot0)
+                out[1] = np.sin(rot0)
+                out[2] = -np.sin(rot1)
+                out[3] = np.cos(rot1)
+        elif corr_mode.literal_value == 2:  # TODO: Is this sensible?
+            def impl(rot0, rot1, out):
+                out[0] = np.cos(rot0)
+                out[1] = np.cos(rot1)
+        elif corr_mode.literal_value == 1:  # TODO: Is this sensible?
+            def impl(rot0, rot1, out):
+                out[0] = np.cos(rot0)
+        else:
+            raise ValueError("Unsupported number of correlations.")
+    else:
+        raise ValueError("Unsupported feed type.")
+
+    return qcjit(impl)
diff --git a/quartical/gains/parallactic_angle/__init__.py b/quartical/gains/parallactic_angle/__init__.py
new file mode 100644
index 00000000..9afc079f
--- /dev/null
+++ b/quartical/gains/parallactic_angle/__init__.py
@@ -0,0 +1,79 @@
+import numpy as np
+from collections import namedtuple
+from quartical.gains.conversion import trig_to_angle
+from quartical.gains.parameterized_gain import ParameterizedGain
+from quartical.gains.parallactic_angle.kernel import (
+    parallactic_angle_params_to_gains
+)
+from quartical.gains.general.flagging import (
+    apply_gain_flags_to_gains,
+    apply_param_flags_to_params
+)
+from quartical.data_handling.angles import _make_parangles
+
+# Overload the default measurement set inputs to include information required
+# to compute the parallactic angles.
+ms_inputs = namedtuple(
+    'ms_inputs', 
+    ParameterizedGain.ms_inputs._fields + \
+        ('RECEPTOR_ANGLE', 'POSITION', 'TIME')
+)
+
+class ParallacticAngle(ParameterizedGain):
+
+    solver = None
+    ms_inputs = ms_inputs
+
+    native_to_converted = (
+        (0, (np.cos,)),
+        (1, (np.sin,))
+    )
+    converted_to_native = (
+        (2, trig_to_angle),
+    )
+    converted_dtype = np.float64
+    native_dtype = np.float64
+
+    def __init__(self, term_name, term_opts):
+
+        super().__init__(term_name, term_opts)
+
+    @classmethod
+    def _make_freq_map(cls, chan_freqs, chan_widths, freq_interval):
+        # Overload gain mapping construction - we evaluate it in every channel.
+        return np.arange(chan_freqs.size, dtype=np.int32)
+
+    @classmethod
+    def make_param_names(cls, correlations):
+        return ["parallactic_angle"]
+
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
+        """Initialise the gains (and parameters)."""
+
+        gains, gain_flags, params, param_flags = super().init_term(
+            term_spec, ref_ant, ms_kwargs, term_kwargs
+        )
+
+        params[...] = _make_parangles(
+            ms_kwargs["TIME"],
+            np.arange(gains.shape[2]),  # TODO: Just for determining n_ant.
+            ms_kwargs["POSITION"],
+            ms_kwargs["RECEPTOR_ANGLE"],
+            meta["FIELD_CENTRE"],
+            "J2000"  # NOTE: Currently hardcoded - will need to be exposed.
+        )[:, None, :, None, :1]  # Assume single value for parallactic angle.
+
+        # Convert the parameters into gains.
+        parallactic_angle_params_to_gains(
+            params,
+            gains,
+            term_kwargs[f"{self.name}_param_freq_map"],
+            feed_type=meta["FEED_TYPE"],
+            corr_mode=gains.shape[-1]
+        )
+
+        # Apply flags to gains and parameters.
+        apply_param_flags_to_params(param_flags, params, 1)
+        apply_gain_flags_to_gains(gain_flags, gains)
+
+        return gains, gain_flags, params, param_flags
diff --git a/quartical/gains/parallactic_angle/kernel.py b/quartical/gains/parallactic_angle/kernel.py
new file mode 100644
index 00000000..177a38ee
--- /dev/null
+++ b/quartical/gains/parallactic_angle/kernel.py
@@ -0,0 +1,71 @@
+# -*- coding: utf-8 -*-
+from numba import njit, literally
+from numba.extending import overload
+from quartical.utils.numba import JIT_OPTIONS
+import quartical.gains.general.factories as factories
+
+
+@njit(**JIT_OPTIONS)
+def parallactic_angle_params_to_gains(
+    params,
+    gains,
+    param_freq_map,
+    feed_type,
+    corr_mode
+):
+    return parallactic_angle_params_to_gains_impl(
+        params,
+        gains,
+        param_freq_map,
+        literally(feed_type),
+        literally(corr_mode)
+    )
+
+
+def parallactic_angle_params_to_gains_impl(
+    params,
+    gains,
+    param_freq_map,
+    feed_type,
+    corr_mode
+):
+    return NotImplementedError
+
+
+@overload(
+    parallactic_angle_params_to_gains_impl,
+    jit_options=JIT_OPTIONS,
+    prefer_literal=True
+)
+def nb_parallactic_angle_params_to_gains_impl(
+    params,
+    gains,
+    param_freq_map,
+    feed_type,
+    corr_mode
+):
+
+    rotmat = factories.rotation_factory(corr_mode, feed_type)
+
+    def impl(
+        params,
+        gains,
+        param_freq_map,
+        feed_type,
+        corr_mode
+    ):
+
+        n_time, n_freq, n_ant, n_dir, n_corr = gains.shape
+
+        for t in range(n_time):
+            for f in range(n_freq):
+                f_m = param_freq_map[f]
+                for a in range(n_ant):
+                    for d in range(n_dir):
+
+                        g = gains[t, f, a, d]
+                        p = params[t, f_m, a, d, 0]
+
+                        rotmat(p, p, g)
+
+    return impl
\ No newline at end of file
diff --git a/quartical/gains/parameterized_gain.py b/quartical/gains/parameterized_gain.py
index 4b50f417..ab68361c 100644
--- a/quartical/gains/parameterized_gain.py
+++ b/quartical/gains/parameterized_gain.py
@@ -174,7 +174,7 @@ def _make_param_freq_chunks(cls, param_freq_map):
     def make_param_names(cls, correlations):
         raise NotImplementedError
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         (_, _, gain_shape, param_shape) = term_spec
diff --git a/quartical/gains/phase/__init__.py b/quartical/gains/phase/__init__.py
index 4b2313ae..857f6323 100644
--- a/quartical/gains/phase/__init__.py
+++ b/quartical/gains/phase/__init__.py
@@ -39,7 +39,7 @@ def make_param_names(cls, correlations):
 
         return [f"phase_{c}" for c in param_corr]
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
diff --git a/quartical/gains/rotation/__init__.py b/quartical/gains/rotation/__init__.py
index 03d03c3d..3ebe5e8e 100644
--- a/quartical/gains/rotation/__init__.py
+++ b/quartical/gains/rotation/__init__.py
@@ -34,7 +34,7 @@ def make_param_names(cls, correlations):
 
         return ["rotation_angle"]
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
diff --git a/quartical/gains/rotation_measure/__init__.py b/quartical/gains/rotation_measure/__init__.py
index 213d8a6b..37f89bb6 100644
--- a/quartical/gains/rotation_measure/__init__.py
+++ b/quartical/gains/rotation_measure/__init__.py
@@ -71,7 +71,7 @@ def make_f_maps(chan_freqs, chan_widths, f_int):
 
         return f_map_arr
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
diff --git a/quartical/gains/tec_and_offset/__init__.py b/quartical/gains/tec_and_offset/__init__.py
index 0360b2c9..6c54b785 100644
--- a/quartical/gains/tec_and_offset/__init__.py
+++ b/quartical/gains/tec_and_offset/__init__.py
@@ -55,7 +55,7 @@ def make_param_names(cls, correlations):
 
         return [n.format(c) for c in param_corr for n in template]
 
-    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs):
+    def init_term(self, term_spec, ref_ant, ms_kwargs, term_kwargs, meta=None):
         """Initialise the gains (and parameters)."""
 
         gains, gain_flags, params, param_flags = super().init_term(
diff --git a/quartical/statistics/logging.py b/quartical/statistics/logging.py
index f0d6021c..8a606593 100644
--- a/quartical/statistics/logging.py
+++ b/quartical/statistics/logging.py
@@ -45,7 +45,9 @@ def log_chisq(pre, post, attrs, block_id=None):
     pre = pre.item()
     post = post.item()
 
-    if np.isfinite(pre) and np.isfinite(post):
+    if pre == 0:  # This may be zero in noise-free tests.
+        fractional_change = 0
+    elif np.isfinite(pre) and np.isfinite(post):
         fractional_change = (post - pre) / pre
     else:
         fractional_change = np.nan
diff --git a/testing/fixtures/data_handling.py b/testing/fixtures/data_handling.py
index 6c6b58a3..0e2ee72e 100644
--- a/testing/fixtures/data_handling.py
+++ b/testing/fixtures/data_handling.py
@@ -2,7 +2,8 @@
 from quartical.data_handling.ms_handler import (read_xds_list,
                                                 preprocess_xds_list)
 from quartical.data_handling.model_handler import add_model_graph
-from quartical.data_handling.angles import make_parangle_xds_list
+from quartical.data_handling.angles import (make_parangle_xds_list,
+                                            assign_parangle_data)
 from quartical.statistics.statistics import make_stats_xds_list
 
 
@@ -23,7 +24,8 @@ def ref_xds_list(read_xds_list_output):
 
 @pytest.fixture(scope="module")
 def preprocessed_xds_list(raw_xds_list, ms_opts):
-    return preprocess_xds_list(raw_xds_list, ms_opts)
+    xdsl = preprocess_xds_list(raw_xds_list, ms_opts)
+    return assign_parangle_data(ms_opts.path, xdsl)
 
 
 @pytest.fixture(scope="module")
diff --git a/testing/tests/gains/test_crosshand_phase_null_v.py b/testing/tests/gains/test_crosshand_phase_null_v.py
new file mode 100644
index 00000000..c92876f4
--- /dev/null
+++ b/testing/tests/gains/test_crosshand_phase_null_v.py
@@ -0,0 +1,172 @@
+from copy import deepcopy
+import pytest
+import numpy as np
+import dask.array as da
+from quartical.calibration.calibrate import add_calibration_graph
+from testing.utils.gains import apply_gains
+
+
+@pytest.fixture(scope="module")
+def opts(base_opts):
+
+    # Don't overwrite base config - instead create a copy and update.
+
+    _opts = deepcopy(base_opts)
+
+    _opts.input_ms.select_corr = [0, 1, 2, 3]
+    _opts.solver.terms = ['G']
+    _opts.solver.iter_recipe = [100]
+    _opts.solver.propagate_flags = False
+    _opts.solver.convergence_criteria = 1e-7
+    _opts.solver.convergence_fraction = 1
+    _opts.solver.threads = 2
+    _opts.G.time_interval = 0
+    _opts.G.type = "crosshand_phase_null_v"
+    _opts.G.solve_per = "array"
+
+    return _opts
+
+
+@pytest.fixture(scope="module")
+def raw_xds_list(read_xds_list_output):
+    # Only use the first xds. This overloads the global fixture.
+    return read_xds_list_output[0][:1]
+
+
+@pytest.fixture(scope="module")
+def true_gain_list(predicted_xds_list):
+
+    gain_list = []
+
+    for xds in predicted_xds_list:
+
+        n_ant = xds.sizes["ant"]
+        utime_chunks = xds.UTIME_CHUNKS
+        n_time = sum(utime_chunks)
+        chan_chunks = xds.chunks["chan"]
+        n_chan = xds.sizes["chan"]
+        n_dir = xds.sizes["dir"]
+        n_corr = xds.sizes["corr"]
+
+        chunking = (len(utime_chunks), chan_chunks, n_ant, n_dir, n_corr)
+
+        bound = np.pi
+
+        da.random.seed(0)
+        amp = da.from_array(np.array([1, 0, 0, 1]))
+        phase = da.random.uniform(size=(1, n_chan, 1, n_dir, n_corr),
+                                  high=bound,
+                                  low=-bound,
+                                  chunks=chunking)
+        phase[..., 1:] = 0  # Only introduce the phase on XX or RR.
+
+        gains = amp[None, None, None, None, :]*da.exp(1j*phase)
+
+        shape = (n_time, n_chan, n_ant, n_dir, n_corr)
+
+        gains = da.broadcast_to(gains[:, :, :1], shape)
+
+        gain_list.append(gains)
+
+    return gain_list
+
+
+@pytest.fixture(scope="module")
+def corrupted_data_xds_list(predicted_xds_list, true_gain_list):
+
+    corrupted_data_xds_list = []
+
+    for xds, gains in zip(predicted_xds_list, true_gain_list):
+
+        n_corr = xds.sizes["corr"]
+
+        ant1 = xds.ANTENNA1.data
+        ant2 = xds.ANTENNA2.data
+        time = xds.TIME.data
+
+        row_inds = \
+            time.map_blocks(lambda x: np.unique(x, return_inverse=True)[1])
+
+        model = da.ones(xds.MODEL_DATA.data.shape, dtype=np.complex128)
+
+        model *= da.from_array([1, 0.1, 0.1, 1])
+
+        data = da.blockwise(apply_gains, ("rfc"),
+                            model, ("rfdc"),
+                            gains, ("rfadc"),
+                            ant1, ("r"),
+                            ant2, ("r"),
+                            row_inds, ("r"),
+                            n_corr, None,
+                            align_arrays=False,
+                            concatenate=True,
+                            dtype=model.dtype)
+
+        corrupted_xds = xds.assign({
+            "DATA": ((xds.DATA.dims), data),
+            "MODEL_DATA": ((xds.MODEL_DATA.dims), model),
+            "FLAG": ((xds.FLAG.dims), da.zeros_like(xds.FLAG.data)),
+            "WEIGHT": ((xds.WEIGHT.dims), da.ones_like(xds.WEIGHT.data))
+            }
+        )
+
+        corrupted_data_xds_list.append(corrupted_xds)
+
+    return corrupted_data_xds_list
+
+
+@pytest.fixture(scope="module")
+def add_calibration_graph_outputs(corrupted_data_xds_list, stats_xds_list,
+                                  solver_opts, chain, output_opts):
+    # Overload this fixture as we need to use the corrupted xdss.
+    return add_calibration_graph(corrupted_data_xds_list, stats_xds_list,
+                                 solver_opts, chain, output_opts)
+
+
+# -----------------------------------------------------------------------------
+
+def test_stokes_v_magnitude(cmp_post_solve_data_xds_list):
+    # Magnitude of stokes V in the corrected data should tend to zero.
+    for xds in cmp_post_solve_data_xds_list:
+        corr_data = xds._CORRECTED_DATA.data
+        res_v = -0.5j*corr_data[...,1] + 0.5j*corr_data[..., 2]
+        np.testing.assert_array_almost_equal(np.abs(res_v), 0)
+
+
+def test_solver_flags(cmp_post_solve_data_xds_list):
+    # The solver should not add addiitonal flags to the test data.
+    for xds in cmp_post_solve_data_xds_list:
+        np.testing.assert_array_equal(xds._FLAG.data, xds.FLAG.data)
+
+
+@pytest.mark.xfail(reason="Ambiguous result - signs may be flipped.")
+def test_gains(cmp_gain_xds_lod, true_gain_list):
+
+    for solved_gain_dict, true_gain in zip(cmp_gain_xds_lod, true_gain_list):
+        solved_gain_xds = solved_gain_dict["G"]
+        solved_gain, solved_flags = da.compute(solved_gain_xds.gains.data,
+                                               solved_gain_xds.gain_flags.data)
+        true_gain = true_gain.compute().copy()  # Make mutable.
+
+        true_gain[np.where(solved_flags)] = 0
+        solved_gain[np.where(solved_flags)] = 0
+
+        # To ensure the missing antenna handling doesn't render this test
+        # useless, check that we have non-zero entries first.
+        assert np.any(solved_gain), "All gains are zero!"
+        np.testing.assert_array_almost_equal(true_gain[:1], solved_gain)
+
+
+def test_gain_flags(cmp_gain_xds_lod):
+
+    for solved_gain_dict in cmp_gain_xds_lod:
+        solved_gain_xds = solved_gain_dict["G"]
+        solved_flags = solved_gain_xds.gain_flags.values
+
+        frows, fchans, fants, fdir = np.where(solved_flags)
+
+        # We know that these antennas are missing in the test data. No other
+        # antennas should have flags.
+        assert set(np.unique(fants)) == {18, 20}
+
+# -----------------------------------------------------------------------------
diff --git a/testing/tests/gains/test_parallactic_angle.py b/testing/tests/gains/test_parallactic_angle.py
new file mode 100644
index 00000000..945f9552
--- /dev/null
+++ b/testing/tests/gains/test_parallactic_angle.py
@@ -0,0 +1,177 @@
+from copy import deepcopy
+import pytest
+import numpy as np
+import dask.array as da
+from quartical.calibration.calibrate import add_calibration_graph
+from testing.utils.gains import apply_gains
+
+
+@pytest.fixture(scope="module")
+def opts(base_opts):
+
+    # Don't overwrite base config - instead create a copy and update.
+
+    _opts = deepcopy(base_opts)
+
+    _opts.input_ms.select_corr = [0, 1, 2, 3]
+    _opts.input_model.apply_p_jones = False
+    _opts.output.apply_p_jones_inv = False
+    _opts.solver.terms = ['G']
+    _opts.solver.iter_recipe = [0]
+    _opts.solver.propagate_flags = False
+    _opts.solver.convergence_criteria = 1e-7
+    _opts.solver.convergence_fraction = 1
+    _opts.solver.threads = 2
+    _opts.G.freq_interval = 0
+    _opts.G.time_interval = 1
+    _opts.G.type = "parallactic_angle"
+
+    return _opts
+
+
+@pytest.fixture(scope="module")
+def raw_xds_list(read_xds_list_output):
+    # Only use the first xds. This overloads the global fixture.
+    return read_xds_list_output[0][:1]
+
+
+@pytest.fixture(scope="module", params=["linear", "circular"])
+def feed_type(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def true_gain_list(predicted_xds_list, parangle_xds_list, feed_type):
+
+    gain_list = []
+
+    for xds, pxds in zip(predicted_xds_list, parangle_xds_list):
+
+        n_ant = xds.sizes["ant"]
+        utime_chunks = xds.UTIME_CHUNKS
+        n_time = sum(utime_chunks)
+        chan_chunks = xds.chunks["chan"]
+        n_chan = xds.sizes["chan"]
+        n_dir = xds.sizes["dir"]
+        n_corr = xds.sizes["corr"]
+
+        shape = (n_time, n_chan, n_ant, n_dir, n_corr)
+        chunking = (utime_chunks, chan_chunks, n_ant, n_dir, n_corr)
+
+        gains = da.zeros(shape, chunks=chunking, dtype=np.complex128)
+        phase = pxds.PARALLACTIC_ANGLES.data[:, None, :, None, :]
+
+        if feed_type == "circular":
+            gains[..., 0] = da.exp(-1j*phase[..., 0])
+            gains[..., 3] = da.exp(1j*phase[..., 1])
+        elif feed_type == "linear":
+            gains[..., 0] = da.cos(phase[..., 0])
+            gains[..., 1] = da.sin(phase[..., 0])
+            gains[..., 2] = -da.sin(phase[..., 1])
+            gains[..., 3] = da.cos(phase[..., 1])
+        else:
+            raise ValueError(f"Unsupported feed type - {feed_type}.")
+
+        gain_list.append(gains)
+
+    return gain_list
+
+
+@pytest.fixture(scope="module")
+def corrupted_data_xds_list(predicted_xds_list, true_gain_list, feed_type):
+
+    corrupted_data_xds_list = []
+
+    for xds, gains in zip(predicted_xds_list, true_gain_list):
+
+        n_corr = xds.sizes["corr"]
+
+        ant1 = xds.ANTENNA1.data
+        ant2 = xds.ANTENNA2.data
+        time = xds.TIME.data
+
+        row_inds = \
+            time.map_blocks(lambda x: np.unique(x, return_inverse=True)[1])
+
+        model = da.ones(xds.MODEL_DATA.data.shape, dtype=np.complex128)
+
+        model *= da.from_array([1, 0.1, 0.1, 1])
+
+        data = da.blockwise(apply_gains, ("rfc"),
+                            model, ("rfdc"),
+                            gains, ("rfadc"),
+                            ant1, ("r"),
+                            ant2, ("r"),
+                            row_inds, ("r"),
+                            n_corr, None,
+                            align_arrays=False,
+                            concatenate=True,
+                            dtype=model.dtype)
+
+        corrupted_xds = xds.assign({
+            "DATA": ((xds.DATA.dims), data),
+            "MODEL_DATA": ((xds.MODEL_DATA.dims), model),
+            "FLAG": ((xds.FLAG.dims), da.zeros_like(xds.FLAG.data)),
+            "WEIGHT": ((xds.WEIGHT.dims), da.ones_like(xds.WEIGHT.data))
+            }
+        )
+
+        corrupted_xds.attrs["FEED_TYPE"] = feed_type  # Change feed type.
+
+        corrupted_data_xds_list.append(corrupted_xds)
+
+    return corrupted_data_xds_list
+
+
+@pytest.fixture(scope="module")
+def add_calibration_graph_outputs(corrupted_data_xds_list, stats_xds_list,
+                                  solver_opts, chain, output_opts):
+    # Overload this fixture as we need to use the corrupted xdss.
+    return add_calibration_graph(corrupted_data_xds_list, stats_xds_list,
+                                 solver_opts, chain, output_opts)
+
+
+# -----------------------------------------------------------------------------
+
+def test_residual_magnitude(cmp_post_solve_data_xds_list):
+    # Magnitude of the residuals should tend to zero.
+    for xds in cmp_post_solve_data_xds_list:
+        np.testing.assert_array_almost_equal(np.abs(xds._RESIDUAL.data), 0)
+
+
+def test_solver_flags(cmp_post_solve_data_xds_list):
+    # The solver should not add addiitonal flags to the test data.
+    for xds in cmp_post_solve_data_xds_list:
+        np.testing.assert_array_equal(xds._FLAG.data, xds.FLAG.data)
+
+
+def test_gains(cmp_gain_xds_lod, true_gain_list):
+
+    for solved_gain_dict, true_gain in zip(cmp_gain_xds_lod, true_gain_list):
+        solved_gain_xds = solved_gain_dict["G"]
+        solved_gain, solved_flags = da.compute(solved_gain_xds.gains.data,
+                                               solved_gain_xds.gain_flags.data)
+        true_gain = true_gain.compute().copy()  # Make mutable.
+
+        true_gain[np.where(solved_flags)] = 0
+        solved_gain[np.where(solved_flags)] = 0
+
+        # To ensure the missing antenna handling doesn't render this test
+        # useless, check that we have non-zero entries first.
+        assert np.any(solved_gain), "All gains are zero!"
+        np.testing.assert_array_almost_equal(true_gain, solved_gain)
+
+
+def test_gain_flags(cmp_gain_xds_lod):
+
+    for solved_gain_dict in cmp_gain_xds_lod:
+        solved_gain_xds = solved_gain_dict["G"]
+        solved_flags = solved_gain_xds.gain_flags.values
+
+        frows, fchans, fants, fdir = np.where(solved_flags)
+
+        # We know that these antennas are missing in the test data. No other
+        # antennas should have flags.
+        assert set(np.unique(fants)) == {18, 20}
+
+# -----------------------------------------------------------------------------