diff --git a/gstlal-calibration/bin/gstlal_compute_strain b/gstlal-calibration/bin/gstlal_compute_strain
index 94824c5d74e2e7dd2e871de461bfef48b03f89ab..b454ec66a8418285fb2b3c6585fa2ca25c095f4d 100755
--- a/gstlal-calibration/bin/gstlal_compute_strain
+++ b/gstlal-calibration/bin/gstlal_compute_strain
@@ -361,6 +361,7 @@ apply_srcq = Config.getboolean("TDCFConfigurations", "applysrcq")
 apply_fs = Config.getboolean("TDCFConfigurations", "applyfs")
 use_coherence = Config.getboolean("TDCFConfigurations", "usecoherence")
 tdcf_default_to_median = Config.getboolean("TDCFConfigurations", "tdcfdefaulttomedian")
+compute_exact_kappas = Config.getboolean("TDCFConfigurations", "computeexactkappas") if "computeexactkappas" in TDCFConfigs else False
 # Boolean for state vector computation
 compute_calib_statevector = Config.getboolean("CalibrationConfigurations", "computecalibstatevector")
 # Boolean for computing factors from filters file
@@ -653,6 +654,62 @@ if factors_from_filters_file or compute_calib_statevector:
 		if factors_from_filters_file and (compute_kappapum or compute_kappauim):
 			raise ValueError("Cannot compute kappa_PUM or kappa_UIM, as the needed EPICS are not contained in the specified filters file.")
 
+	try:
+		EP25_real = float(filters["EP25_real"])
+		EP25_imag = float(filters["EP25_imag"])
+		EP26_real = float(filters["EP26_real"])
+		EP26_imag = float(filters["EP26_imag"])
+		EP27_real = float(filters["EP27_real"])
+		EP27_imag = float(filters["EP27_imag"])
+		EP28_real = float(filters["EP28_real"])
+		EP28_imag = float(filters["EP28_imag"])
+		EP29_real = float(filters["EP29_real"])
+		EP29_imag = float(filters["EP29_imag"])
+		EP30_real = float(filters["EP30_real"])
+		EP30_imag = float(filters["EP30_imag"])
+		EP31_real = float(filters["EP31_real"])
+		EP31_imag = float(filters["EP31_imag"])
+		EP32_real = float(filters["EP32_real"])
+		EP32_imag = float(filters["EP32_imag"])
+		EP33_real = float(filters["EP33_real"])
+		EP33_imag = float(filters["EP33_imag"])
+		EP34_real = float(filters["EP34_real"])
+		EP34_imag = float(filters["EP34_imag"])
+		EP35_real = float(filters["EP35_real"])
+		EP35_imag = float(filters["EP35_imag"])
+		EP36_real = float(filters["EP36_real"])
+		EP36_imag = float(filters["EP36_imag"])
+		EP37_real = float(filters["EP37_real"])
+		EP37_imag = float(filters["EP37_imag"])
+		EP38_real = float(filters["EP38_real"])
+		EP38_imag = float(filters["EP38_imag"])
+		EP39_real = float(filters["EP39_real"])
+		EP39_imag = float(filters["EP39_imag"])
+		EP40_real = float(filters["EP40_real"])
+		EP40_imag = float(filters["EP40_imag"])
+		EP41_real = float(filters["EP41_real"])
+		EP41_imag = float(filters["EP41_imag"])
+		EP42_real = float(filters["EP42_real"])
+		EP42_imag = float(filters["EP42_imag"])
+	except:
+		if compute_exact_kappas and (compute_kappatst or compute_kappapum or compute_kappauim or compute_kappac or compute_fcc or compute_srcq or compute_fs):
+			raise ValueError("Cannot use ComputeExactKappas, as the needed EPICS (25 - 42) are not in the filters file")
+
+if compute_exact_kappas and (compute_kappatst or compute_kappapum or compute_kappauim or compute_kappac or compute_fcc or compute_srcq or compute_fs):
+	try:
+		EP43_real = float(filters["EP43_real"])
+		EP43_imag = float(filters["EP43_imag"])
+		EP44_real = float(filters["EP44_real"])
+		EP44_imag = float(filters["EP44_imag"])
+		EP45_real = float(filters["EP45_real"])
+		EP45_imag = float(filters["EP45_imag"])
+		EP46_real = float(filters["EP46_real"])
+		EP46_imag = float(filters["EP46_imag"])
+		EP47_real = float(filters["EP47_real"])
+		EP47_imag = float(filters["EP47_imag"])
+	except:
+		raise ValueError("Cannot use ComputeExactKappas, as EPICS 43 - 47 are not in the filters file.")
+
 # Load the high-pass filters for h(t)
 if "ctrl_highpass" in filters:
 	act_highpass = filters["ctrl_highpass"]
@@ -987,6 +1044,8 @@ if (ChannelNames["linewitnesschannellist"] != "None") if "linewitnesschannellist
 			line_witness_frequency_channel_list.append([])
 			for j in range(0, len(line_witness_freqs[i])):
 				line_witness_frequency_channel_list[i].append(None)
+else:
+	line_witness_frequency_channel_list = []
 
 if (ChannelNames["powerlineschannel"] not in headkeys) if remove_power_lines else False:
 	channel_list.append((instrument, ChannelNames["powerlineschannel"]))
@@ -1163,7 +1222,7 @@ if use_coherence:
 		else:
 			pcaly_line4_coh = pcaly_line1_coh
 
-if compute_kappatst or compute_kappapu or compute_kappac or compute_fcc or compute_fs or compute_srcq or (remove_cal_lines and "tstexc_linefreq" in act_line_removal_dict.keys()):
+if compute_kappatst or compute_exact_kappas or compute_kappapu or compute_kappac or compute_fcc or compute_fs or compute_srcq or (remove_cal_lines and "tstexc_linefreq" in act_line_removal_dict.keys()):
 	tstexccaps = "audio/x-raw, format=F64LE, rate=%d" % tst_exc_sr
 	tstexc = calibration_parts.caps_and_progress(pipeline, head_dict["tstexc"], tstexccaps, "tstexc")
 	tstexc = pipeparts.mktee(pipeline, tstexc)
@@ -1174,7 +1233,7 @@ if compute_kappatst or compute_kappapu or compute_kappac or compute_fcc or compu
 				line_witness_channel_list[i][j] = tstexc
 	if "tstexc_linefreq" in act_line_removal_dict:
 		act_line_removal_dict["tstexc_linefreq"][0] = tstexc
-if compute_kappapum or (not compute_kappapu and (compute_kappac or compute_fcc or compute_fs or compute_srcq)) or (remove_cal_lines and "pumexc_linefreq" in act_line_removal_dict.keys()):
+if compute_kappapum or compute_exact_kappas or (not compute_kappapu and (compute_kappac or compute_fcc or compute_fs or compute_srcq)) or (remove_cal_lines and "pumexc_linefreq" in act_line_removal_dict.keys()):
 	pumexccaps = "audio/x-raw, format=F64LE, rate=%d" % pum_exc_sr
 	pumexc = calibration_parts.caps_and_progress(pipeline, head_dict["pumexc"], pumexccaps, "pumexc")
 	pumexc = pipeparts.mktee(pipeline, pumexc)
@@ -1185,7 +1244,7 @@ if compute_kappapum or (not compute_kappapu and (compute_kappac or compute_fcc o
 				line_witness_channel_list[i][j] = pumexc
 	if "pumexc_linefreq" in act_line_removal_dict:
 		act_line_removal_dict["pumexc_linefreq"][0] = pumexc
-if compute_kappauim or (not compute_kappapu and (compute_kappac or compute_fcc or compute_fs or compute_srcq)) or (remove_cal_lines and "uimexc_linefreq" in act_line_removal_dict.keys()):
+if compute_kappauim or compute_exact_kappas or (not compute_kappapu and (compute_kappac or compute_fcc or compute_fs or compute_srcq)) or (remove_cal_lines and "uimexc_linefreq" in act_line_removal_dict.keys()):
 	uimexccaps = "audio/x-raw, format=F64LE, rate=%d" % uim_exc_sr
 	uimexc = calibration_parts.caps_and_progress(pipeline, head_dict["uimexc"], uimexccaps, "uimexc")
 	uimexc = pipeparts.mktee(pipeline, uimexc)
@@ -1201,9 +1260,189 @@ if compute_kappapu:
 	darmexc = calibration_parts.caps_and_progress(pipeline, head_dict["darmexc"], hoft_caps, "darmexc")
 
 # Set up computations for kappa_tst, kappa_pum, kappa_uim, kappa_pu,kappa_c, f_cc, f_s, and Q, if applicable
-if compute_kappatst or compute_kappapum or compute_kappauim or compute_kappapu or compute_kappac or compute_fcc or compute_srcq or compute_fs:
+if compute_exact_kappas and (compute_kappatst or compute_kappapum or compute_kappauim or compute_kappac or compute_fcc or compute_srcq or compute_fs):
+	# Use the "exact" algebraic solution of P1900052, Section 5.2.6.
+	# Do all the demodulation up front.  We need DARM_ERR and either Pcal or an actuator
+	# injection channel at each calibration line frequency used for the TDCFs.  We will
+	# also take the ratios injection(f) / DARM_ERR(f) at all the calibration lines and
+	# apply the gating and smoothing here, before computing the TDCFs.
+
+	# First, the actuation Pcal line, starting with DARM_ERR.
+	derr_at_act_pcal_freq = calibration_parts.demodulate(pipeline, derrtee, act_pcal_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["pcal1_linefreq"] if "pcal1_linefreq" in head_dict else None)
+	# Now Pcal
+	pcal_at_act_pcal_freq = calibration_parts.demodulate(pipeline, pcaltee, act_pcal_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, prefactor_real = pcal_sign * pcal_corr_at_act_freq_real, prefactor_imag = pcal_sign * pcal_corr_at_act_freq_imag, freq_update = [head_dict["pcal1_linefreq"], head_dict["pcal1_line_corr_real"], head_dict["pcal1_line_corr_imag"]] if "pcal1_linefreq" in head_dict else None)
+	pcal_at_act_pcal_freq = pipeparts.mktee(pipeline, pcal_at_act_pcal_freq)
+	if remove_cal_lines:
+		# This will save having to demodulate it again
+		pcal_line_removal_dict["pcal1_linefreq"][0] = pcal_at_act_pcal_freq
+		pcal_line_removal_dict["pcal1_linefreq"][4] = True
+
+	# Take the ratio X1 = Pcal(f1) / DARM_ERR(f1)
+	X1 = calibration_parts.complex_division(pipeline, pcal_at_act_pcal_freq, derr_at_act_pcal_freq)
+
+	# Now the gating and smoothing
+	if use_coherence:
+		X1 = calibration_parts.mkgate(pipeline, X1, pcaly_line1_coh, coherence_unc_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True, name = 'pcal_line1_gate')
+		X1 = calibration_parts.smooth_complex_kappas(pipeline, X1, EP43_real, EP43_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+	else:
+		var = 0.2 * abs(EP43_real + 1j * EP43_imag)
+		X1 = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, X1, var, var, EP43_real, EP43_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+
+	# The sensing Pcal line (~400 Hz), starting with DARM_ERR.
+	derr_at_opt_gain_freq = calibration_parts.demodulate(pipeline, derrtee, opt_gain_fcc_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["pcal2_linefreq"] if "pcal2_linefreq" in head_dict else None)
+	# Now Pcal
+	pcal_at_opt_gain_freq = calibration_parts.demodulate(pipeline, pcaltee, opt_gain_fcc_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, prefactor_real = pcal_sign * pcal_corr_at_opt_gain_fcc_freq_real, prefactor_imag = pcal_sign * pcal_corr_at_opt_gain_fcc_freq_imag, freq_update = [head_dict["pcal2_linefreq"], head_dict["pcal2_line_corr_real"], head_dict["pcal2_line_corr_imag"]] if "pcal2_linefreq" in head_dict else None)
+	if remove_cal_lines:
+		# This will save having to demodulate it again
+		pcal_at_opt_gain_freq = pipeparts.mktee(pipeline, pcal_at_opt_gain_freq)
+		pcal_line_removal_dict["pcal2_linefreq"][0] = pcal_at_opt_gain_freq
+		pcal_line_removal_dict["pcal2_linefreq"][4] = True
+
+	# Take the ratio X2 = Pcal(f2) / DARM_ERR(f2)
+	X2 = calibration_parts.complex_division(pipeline, pcal_at_opt_gain_freq, derr_at_opt_gain_freq)
+
+	# Now the gating and smoothing
+	if use_coherence:
+		X2 = calibration_parts.mkgate(pipeline, X2, pcaly_line2_coh, coherence_unc_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True, name = 'pcal_line2_gate')
+		X2 = calibration_parts.smooth_complex_kappas(pipeline, X2, EP44_real, EP44_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+	else:
+		var = 0.2 * abs(EP44_real + 1j * EP44_imag)
+		X2 = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, X2, var, var, EP44_real, EP44_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+
+	# Now the ESD/TST/L3 actuator line, starting with DARM_ERR.
+	derr_at_esd_act_freq = calibration_parts.demodulate(pipeline, derrtee, esd_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["tstexc_linefreq"] if "tstexc_linefreq" in head_dict else None)
+	# Now the TST excitation channel
+	tstexc_at_esd_act_freq = calibration_parts.demodulate(pipeline, tstexc, esd_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["tstexc_linefreq"] if "tstexc_linefreq" in head_dict else None)
+	if "tstexc_linefreq" in act_line_removal_dict.keys():
+		tstexc_at_esd_act_freq = pipeparts.mktee(pipeline, tstexc_at_esd_act_freq)
+		act_line_removal_dict["tstexc_linefreq"][0] = tstexc_at_esd_act_freq
+		act_line_removal_dict["tstexc_linefreq"][4] = True
+
+	# Take the ratio X_T = Pcal(f_T) / DARM_ERR(f_T)
+	X_T = calibration_parts.complex_division(pipeline, tstexc_at_esd_act_freq, derr_at_esd_act_freq)
+
+	# Now the gating and smoothing
+	if use_coherence:
+		X_T = calibration_parts.mkgate(pipeline, X_T, sus_line3_coh, coherence_unc_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True, name = 'sus_line3_gate')
+		X_T = calibration_parts.smooth_complex_kappas(pipeline, X_T, EP45_real, EP45_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+	else:
+		var = 0.2 * abs(EP45_real + 1j * EP45_imag)
+		X_T = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, X_T, var, var, EP45_real, EP45_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+
+	# Now the PUM/L2 actuator line, starting with DARM_ERR.
+	derr_at_pum_act_freq = calibration_parts.demodulate(pipeline, derrtee, pum_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["pumexc_linefreq"] if "pumexc_linefreq" in head_dict else None)
+	# Now the PUM excitation channel
+	pumexc_at_pum_act_freq = calibration_parts.demodulate(pipeline, pumexc, pum_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["pumexc_linefreq"] if "pumexc_linefreq" in head_dict else None)
+	if "pumexc_linefreq" in act_line_removal_dict.keys():
+		pumexc_at_pum_act_freq = pipeparts.mktee(pipeline, pumexc_at_pum_act_freq)
+		act_line_removal_dict["pumexc_linefreq"][0] = pumexc_at_pum_act_freq
+		act_line_removal_dict["pumexc_linefreq"][4] = True
+
+	# Take the ratio X_P = Pcal(f_P) / DARM_ERR(f_P)
+	X_P = calibration_parts.complex_division(pipeline, pumexc_at_pum_act_freq, derr_at_pum_act_freq)
+
+	# Now the gating and smoothing
+	if use_coherence:
+		X_P = calibration_parts.mkgate(pipeline, X_P, sus_line2_coh, coherence_unc_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True, name = 'sus_line2_gate')
+		X_P = calibration_parts.smooth_complex_kappas(pipeline, X_P, EP46_real, EP46_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+	else:
+		var = 0.2 * abs(EP46_real + 1j * EP46_imag)
+		X_P = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, X_P, var, var, EP46_real, EP46_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+
+	# Now the UIM/L3 actuator line, starting with DARM_ERR.
+	derr_at_uim_act_freq = calibration_parts.demodulate(pipeline, derrtee, uim_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["uimexc_linefreq"] if "uimexc_linefreq" in head_dict else None)
+	# Now the UIM excitation channel
+	uimexc_at_uim_act_freq = calibration_parts.demodulate(pipeline, uimexc, uim_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["uimexc_linefreq"] if "uimexc_linefreq" in head_dict else None)
+	if "uimexc_linefreq" in act_line_removal_dict.keys():
+		uimexc_at_uim_act_freq = pipeparts.mktee(pipeline, uimexc_at_uim_act_freq)
+		act_line_removal_dict["uimexc_linefreq"][0] = uimexc_at_uim_act_freq
+		act_line_removal_dict["uimexc_linefreq"][4] = True
+
+	# Take the ratio X_U = Pcal(f_U) / DARM_ERR(f_U)
+	X_U = calibration_parts.complex_division(pipeline, uimexc_at_uim_act_freq, derr_at_uim_act_freq)
+
+	# Now the gating and smoothing
+	if use_coherence:
+		X_U = calibration_parts.mkgate(pipeline, X_U, sus_line1_coh, coherence_unc_threshold, attack_length = kappa_gate_attack_length, hold_length = kappa_gate_hold_length, invert_control = True, name = 'sus_line1_gate')
+		X_U = calibration_parts.smooth_complex_kappas(pipeline, X_U, EP47_real, EP47_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+	else:
+		var = 0.2 * abs(EP47_real + 1j * EP47_imag)
+		X_U = calibration_parts.smooth_complex_kappas_no_coherence(pipeline, X_U, var, var, EP47_real, EP47_imag, median_smoothing_samples, factors_average_samples, tdcf_default_to_median, filter_latency_factor)
+
+	# Finally, compute the kappas
+	if factors_from_filters_file:
+		X_list = [X1, X2, X_T, X_P, X_U]
+		freq_list = [act_pcal_line_freq, opt_gain_fcc_line_freq, esd_act_line_freq, pum_act_line_freq, uim_act_line_freq]
+		EPICS_list = [EP11_real, EP11_imag, EP25_real, EP25_imag, EP26_real, EP26_imag, EP27_real, EP27_imag, EP6_real, EP6_imag, EP28_real, EP28_imag, EP29_real, EP29_imag, EP30_real, EP30_imag, EP31_real, EP31_imag, EP32_real, EP32_imag, EP33_real, EP33_imag, EP34_real, EP34_imag, EP35_real, EP35_imag, EP36_real, EP36_imag, EP37_real, EP37_imag, EP38_real, EP38_imag, EP39_real, EP39_imag, EP40_real, EP40_imag, EP41_real, EP41_imag, EP42_real, EP42_imag]
+
+		[ktst, kpum, kuim, tau_tst, tau_pum, tau_uim, kc, fcc, fs_squared, fs_over_Q] = calibration_parts.compute_exact_kappas_from_filters_file(pipeline, X_list, freq_list, EPICS_list)
+
+	#else: use non-existent channels in the raw frames
+
+	# A little bit of processing is necessary to ensure backwards compatibility
+	if compute_kappatst:
+		smooth_ktsttee = pipeparts.mktee(pipeline, calibration_parts.mkmultiplier(pipeline, [pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, ktst, matrix = [[1.0, 0.0]])), pipeparts.mkgeneric(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, tau_tst, matrix = [[0.0, 2.0 * numpy.pi * esd_act_line_freq]])), "cexp")]))
+		smooth_ktstR, smooth_ktstI = calibration_parts.split_into_real(pipeline, smooth_ktsttee)
+		smooth_ktstRtee = pipeparts.mktee(pipeline, smooth_ktstR)
+		smooth_ktstItee = pipeparts.mktee(pipeline, smooth_ktstI)
+		smooth_ktstR_nogate = smooth_ktstRtee
+		smooth_ktstI_nogate = smooth_ktstItee
+	else:
+		pipeparts.mkfakesink(pipeline, ktst)
+		pipeparts.mkfakesink(pipeline, tau_tst)
+
+	if compute_kappapum:
+		smooth_kpumtee = pipeparts.mktee(pipeline, calibration_parts.mkmultiplier(pipeline, [pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, kpum, matrix = [[1.0, 0.0]])), pipeparts.mkgeneric(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, tau_pum, matrix = [[0.0, 2.0 * numpy.pi * pum_act_line_freq]])), "cexp")]))
+		smooth_kpumR, smooth_kpumI = calibration_parts.split_into_real(pipeline, smooth_kpumtee)
+		smooth_kpumRtee = pipeparts.mktee(pipeline, smooth_kpumR)
+		smooth_kpumItee = pipeparts.mktee(pipeline, smooth_kpumI)
+		smooth_kpumR_nogate = smooth_kpumRtee
+		smooth_kpumI_nogate = smooth_kpumItee
+	else:
+		pipeparts.mkfakesink(pipeline, kpum)
+		pipeparts.mkfakesink(pipeline, tau_pum)
+
+	if compute_kappauim:
+		smooth_kuimtee = pipeparts.mktee(pipeline, calibration_parts.mkmultiplier(pipeline, [pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, kuim, matrix = [[1.0, 0.0]])), pipeparts.mkgeneric(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, tau_uim, matrix = [[0.0, 2.0 * numpy.pi * uim_act_line_freq]])), "cexp")]))
+		smooth_kuimR, smooth_kuimI = calibration_parts.split_into_real(pipeline, smooth_kuimtee)
+		smooth_kuimRtee = pipeparts.mktee(pipeline, smooth_kuimR)
+		smooth_kuimItee = pipeparts.mktee(pipeline, smooth_kuimI)
+		smooth_kuimR_nogate = smooth_kuimRtee
+		smooth_kuimI_nogate = smooth_kuimItee
+	else:
+		pipeparts.mkfakesink(pipeline, kuim)
+		pipeparts.mkfakesink(pipeline, tau_uim)
+
+	if compute_kappac:
+		smooth_kctee = smooth_kc_nogate = kc
+	else:
+		pipeparts.mkfakesink(pipeline, kc)
+
+	if compute_fcc:
+		smooth_fcctee = smooth_fcc_nogate = fcc
+	else:
+		pipeparts.mkfakesink(pipeline, fcc)
+
+	if compute_fs:
+		smooth_fs_squared_nogate = smooth_fs_squared = fs_squared
+	else:
+		pipeparts.mkfakesink(pipeline, fs_squared)
+
+	if compute_srcq:
+		# fs / Q is a real-valued quantity, but fs and Q are not necessarily, so we compute a complex-valued 1/Q.
+		smooth_srcQ_inv = smooth_srcQ_inv_nogate = pipeparts.mktee(pipeline, calibration_parts.mkmultiplier(pipeline, [pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, fs_over_Q, matrix = [[1.0, 0.0]])), calibration_parts.mkpow(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, fs_squared, matrix = [[1.0, 0.0]])), exponent = -0.5)]))
+		# We need a real-valued version of Q^(-1) to write the the frames.  If we have
+		# an optical spring, the computed Q^(-1) is imaginary, and if we have an optical
+		# antispring, the computed Q^(-1) is real.  To get a sensible real value either
+		# way AND INCLUDE ANY MINUS SIGN, we use Re(Q^(-1)) + Im(Q^(-1)).
+		smooth_srcQ_inv_real = smooth_srcQ_inv_nogate = pipeparts.mktee(pipeline, calibration_parts.mkadder(pipeline, calibration_parts.list_srcs(pipeline, pipeparts.mkgeneric(pipeline, smooth_srcQ_inv, "creal"), pipeparts.mkgeneric(pipeline, smooth_srcQ_inv, "cimag"))))
+
+	else:
+		pipeparts.mkfakesink(pipeline, fs_over_Q)
 
-	# demodulate the PCAL channel and apply the PCAL correction factor at the DARM actuation line frequency
+elif compute_kappatst or compute_kappapum or compute_kappauim or compute_kappapu or compute_kappac or compute_fcc or compute_srcq or compute_fs:
+	# Use the solution in T1700106.
+	# demodulate the PCAL channel and apply the PCAL correction factor at the actuation Pcal line frequency
 	pcal_at_act_pcal_freq = calibration_parts.demodulate(pipeline, pcaltee, act_pcal_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, prefactor_real = pcal_sign * pcal_corr_at_act_freq_real, prefactor_imag = pcal_sign * pcal_corr_at_act_freq_imag, freq_update = [head_dict["pcal1_linefreq"], head_dict["pcal1_line_corr_real"], head_dict["pcal1_line_corr_imag"]] if "pcal1_linefreq" in head_dict else None)
 	pcal_at_act_pcal_freq = pipeparts.mktee(pipeline, pcal_at_act_pcal_freq)
 	if remove_cal_lines:
@@ -1211,7 +1450,7 @@ if compute_kappatst or compute_kappapum or compute_kappauim or compute_kappapu o
 		pcal_line_removal_dict["pcal1_linefreq"][0] = pcal_at_act_pcal_freq
 		pcal_line_removal_dict["pcal1_linefreq"][4] = True
 
-	# demodulate DARM_ERR at the ~30 Hz pcal line frequency
+	# demodulate DARM_ERR at the actuation pcal line frequency
 	derr_at_act_pcal_freq = calibration_parts.demodulate(pipeline, derrtee, act_pcal_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["pcal1_linefreq"] if "pcal1_linefreq" in head_dict else None)
 	if dewhitening:
 		# dewhiten DARM_ERR at the ~30 Hz pcal line frequency
@@ -1274,9 +1513,10 @@ if compute_kappatst or compute_kappapum or compute_kappauim or compute_kappapu o
 if compute_kappapum or (not compute_kappapu and (compute_kappac or compute_fcc or compute_fs or compute_srcq)):
 	# demodulate the PUM excitation channel at the PUM actuation line frequency
 	pumexc_at_pum_act_freq = calibration_parts.demodulate(pipeline, pumexc, pum_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["pumexc_linefreq"] if "pumexc_linefreq" in head_dict else None)
-	pumexc_at_pum_act_freq = pipeparts.mktee(pipeline, pumexc_at_pum_act_freq)
-	act_line_removal_dict["pumexc_linefreq"][0] = pumexc_at_pum_act_freq
-	act_line_removal_dict["pumexc_linefreq"][4] = True
+	if "pumexc_linefreq" in act_line_removal_dict.keys():
+		pumexc_at_pum_act_freq = pipeparts.mktee(pipeline, pumexc_at_pum_act_freq)
+		act_line_removal_dict["pumexc_linefreq"][0] = pumexc_at_pum_act_freq
+		act_line_removal_dict["pumexc_linefreq"][4] = True
 
 	# demodulate DARM_ERR at the PUM actuation line frequency
 	derr_at_pum_act_freq = calibration_parts.demodulate(pipeline, derrtee, pum_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["pumexc_linefreq"] if "pumexc_linefreq" in head_dict else None)
@@ -1325,9 +1565,10 @@ if compute_kappauim or (not compute_kappapu and (compute_kappac or compute_fcc o
 	# Demodulate DARM_ERR and the UIM excitation channel at the UIM actuation line frequency
 	derr_at_uim_act_freq = calibration_parts.demodulate(pipeline, derrtee, uim_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["uimexc_linefreq"] if "uimexc_linefreq" in head_dict else None)
 	uimexc_at_uim_act_freq = calibration_parts.demodulate(pipeline, uimexc, uim_act_line_freq, td, compute_factors_sr, demodulation_filter_time, filter_latency_factor, freq_update = head_dict["uimexc_linefreq"] if "uimexc_linefreq" in head_dict else None)
-	uimexc_at_uim_act_freq = pipeparts.mktee(pipeline, uimexc_at_uim_act_freq)
-	act_line_removal_dict["uimexc_linefreq"][0] = uimexc_at_uim_act_freq
-	act_line_removal_dict["uimexc_linefreq"][4] = True
+	if "uimexc_linefreq" in act_line_removal_dict.keys():
+		uimexc_at_uim_act_freq = pipeparts.mktee(pipeline, uimexc_at_uim_act_freq)
+		act_line_removal_dict["uimexc_linefreq"][0] = uimexc_at_uim_act_freq
+		act_line_removal_dict["uimexc_linefreq"][4] = True
 
 	# Compute kappa_uim, either using reference factors from the filters file or reading them from EPICS channels
 	if not factors_from_filters_file:
diff --git a/gstlal-calibration/python/calibration_parts.py b/gstlal-calibration/python/calibration_parts.py
index a4b48f0b22eac83ed5767e00d216b631ef26ff07..ab44ee3fa771a7a50d25783d00b2a6488ef9881d 100644
--- a/gstlal-calibration/python/calibration_parts.py
+++ b/gstlal-calibration/python/calibration_parts.py
@@ -958,6 +958,193 @@ def compute_kappaa(pipeline, afctrl, EP4, EP5):
 
 	return ka
 
+def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS):
+
+	#
+	# See P1900052, Section 5.2.6 for details.  All constants are contained in the list
+	# variable EPICS.  The variable freqs is a list containing calibration line
+	# frequencies, stored in the order f1, f2, fT, fP, fU, i.e., the Pcal lines come
+	# first, and then the actuator lines.  All other quantities evaluated at the
+	# calibration lines are stored in the same order.  The list variable X contains the
+	# ratios X[i] = injection(f_i) / d_err(f_i) for each calibration line frequency.
+	#
+
+	print("starting exact kappas")
+	kappas = []
+	num_lines = len(freqs)
+	num_stages = num_lines - 2 # Stages of actuation (currently 3)
+	MV_matrix = list(numpy.zeros(2 * num_stages * (2 * num_stages + 1)))
+	Yreal = []
+	Yimag = []
+	CAX = []
+	CAXreal = []
+	CAXimag = []
+	Gresreal = []
+	Gresimag = []
+	kappas = []
+
+	for i in range(num_lines):
+		if i < 2:
+			# Then it's a Pcal line
+			Y = pipeparts.mktee(pipeline, complex_audioamplify(pipeline, X[i], EPICS[2 * (1 + num_stages) * i], EPICS[2 * (1 + num_stages) * i + 1]))
+			Yreal.append(pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkgeneric(pipeline, Y, "creal"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0", name = "capsfilter_Yreal_%d" % i)))
+			Yimag.append(pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkgeneric(pipeline, Y, "cimag"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0", name = "capsfilter_Yimag_%d" % i)))
+		else:
+			# It's an actuator line
+			CAX.append(pipeparts.mktee(pipeline, complex_audioamplify(pipeline, X[i], EPICS[2 * (1 + num_stages) * i], EPICS[2 * (1 + num_stages) * i + 1])))
+			CAXreal.append(pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkgeneric(pipeline, CAX[-1], "creal"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0", name = "capsfilter_CAXreal_%d" % i)))
+			CAXimag.append(pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkgeneric(pipeline, CAX[-1], "cimag"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0", name = "capsfilter_CAXimag_%d" % i)))
+
+	# Let's start by computing the V's of Eqs. 5.2.78 and 5.2.79
+	for j in range(num_stages):
+		factor1 = pow(freqs[0], -2) - pow(freqs[2 + j], -2)
+		factor2 = pow(freqs[2 + j], -2) - pow(freqs[1], -2)
+		factor3 = pow(freqs[0], 2) - pow(freqs[2 + j], 2)
+		factor4 = freqs[0] * (freqs[2 + j] - freqs[1])
+		Vj = mkadder(pipeline, list_srcs(pipeline, pipeparts.mkaudioamplify(pipeline, Yreal[1], factor1), pipeparts.mkaudioamplify(pipeline, Yreal[0], factor2)))
+		Vjplus3 = mkadder(pipeline, list_srcs(pipeline, pipeparts.mkaudioamplify(pipeline, Yimag[1], factor3), pipeparts.mkaudioamplify(pipeline, Yimag[0], factor4)))
+		MV_matrix[j] = Vj
+		MV_matrix[num_stages + j] = Vjplus3
+
+	# Now let's compute the elements of the matrix M, given by Eqs. 5.2.70 - 5.2.77
+	# Many of the elements are constant, so make a stream of ones to multiply
+	ones = pipeparts.mktee(pipeline, mkpow(pipeline, Yreal[0], exponent = 0.0))
+	print("exact kappas 10")
+	for j in range(num_stages):
+		# Time-dependent matrix elements
+		factor = pow(freqs[0], -2) - pow(freqs[1], -2)
+		addend = (pow(freqs[0], -2) - pow(freqs[2 + j], -2)) * EPICS[2 * ((1 + num_stages) + 1 + j)] + (pow(freqs[2 + j], -2) - pow(freqs[1], -2)) * EPICS[2 * (1 + j)] - (pow(freqs[0], -2) - pow(freqs[1], -2)) * EPICS[2 * ((2 + j) * (1 + num_stages) + 1 + j)]
+		Mjj = pipeparts.mkgeneric(pipeline, pipeparts.mkaudioamplify(pipeline, CAXreal[j], factor), "lal_add_constant", value = addend)
+		factor = -2.0 * numpy.pi * freqs[2 + j] * (pow(freqs[0], -2) - pow(freqs[1], -2))
+		addend = -2.0 * numpy.pi * freqs[1] * (pow(freqs[0], -2) - pow(freqs[2 + j], -2)) * EPICS[1 + 2 * ((1 + num_stages) + 1 + j)] - 2.0 * numpy.pi * freqs[0] * (pow(freqs[2 + j], -2) - pow(freqs[1], -2)) * EPICS[1 + 2 * (1 + j)] + 2.0 * numpy.pi * freqs[2 + j] * (pow(freqs[0], -2) - pow(freqs[1], -2)) * EPICS[1 + 2 * ((2 + j) * (1 + num_stages) + 1 + j)]
+		Mjjplus3 = pipeparts.mkgeneric(pipeline, pipeparts.mkaudioamplify(pipeline, CAXimag[j], factor), "lal_add_constant", value = addend)
+		factor = -freqs[2 + j] * (pow(freqs[1], 2) - pow(freqs[0], 2))
+		addend = freqs[1] * (pow(freqs[0], 2) - pow(freqs[2 + j], 2)) * EPICS[1 + 2 * ((1 + num_stages) + 1 + j)] + freqs[0] * (pow(freqs[2 + j], 2) - pow(freqs[1], 2)) * EPICS[1 + 2 * (1 + j)] + freqs[2 + j] * (pow(freqs[1], 2) - pow(freqs[0], 2)) * EPICS[1 + 2 * ((2 + j) * (1 + num_stages) + 1 + j)]
+		Mjplus3j = pipeparts.mkgeneric(pipeline, pipeparts.mkaudioamplify(pipeline, CAXimag[j], factor), "lal_add_constant", value = addend)
+		factor = -2.0 * numpy.pi * pow(freqs[2 + j], 2) * (pow(freqs[1], 2) - pow(freqs[0], 2))
+		addend = 2.0 * numpy.pi * pow(freqs[1], 2) * (pow(freqs[0], 2) - pow(freqs[2 + j], 2)) * EPICS[2 * ((1 + num_stages) + 1 + j)] + 2.0 * numpy.pi * pow(freqs[0], 2) * (pow(freqs[2 + j], 2) - pow(freqs[1], 2)) * EPICS[2 * (1 + j)] + 2.0 * numpy.pi * pow(freqs[2 + j], 2) * (pow(freqs[1], 2) - pow(freqs[0], 2)) * EPICS[2 * ((2 + j) * (1 + num_stages) + 1 + j)]
+		Mjplus3jplus3 = pipeparts.mkgeneric(pipeline, pipeparts.mkaudioamplify(pipeline, CAXreal[j], factor), "lal_add_constant", value = addend)
+		# Add these into the matrix
+		MV_matrix[(1 + j) * 2 * num_stages + j] = Mjj
+		MV_matrix[(1 + j) * 2 * num_stages + num_stages + j] = Mjjplus3
+		MV_matrix[(1 + num_stages + j) * 2 * num_stages + j] = Mjplus3j
+		MV_matrix[(1 + num_stages + j) * 2 * num_stages + num_stages + j] = Mjplus3jplus3
+
+		# Constant matrix elements
+		knotequalj = range(num_stages)
+		knotequalj.remove(j)
+		for k in knotequalj:
+			factor = (pow(freqs[0], -2) - pow(freqs[2 + j], -2)) * EPICS[2 * ((1 + num_stages) + 1 + k)] + (pow(freqs[2 + j], -2) - pow(freqs[1], -2)) * EPICS[2 * (1 + k)] - (pow(freqs[0], -2) - pow(freqs[1], -2)) * EPICS[2 * ((2 + j) * (1 + num_stages) + 1 + k)]
+			Mjk = pipeparts.mkaudioamplify(pipeline, ones, factor)
+			factor = -2.0 * numpy.pi * freqs[1] * (pow(freqs[0], -2) - pow(freqs[2 + j], -2)) * EPICS[1 + 2 * ((1 + num_stages) + 1 + k)] - 2.0 * numpy.pi * freqs[0] * (pow(freqs[2 + j], -2) - pow(freqs[1], -2)) * EPICS[1 + 2 * (1 + k)] + 2.0 * numpy.pi * freqs[2 + j] * (pow(freqs[0], -2) - pow(freqs[1], -2)) * EPICS[1 + 2 * ((2 + j) * (1 + num_stages) + 1 + k)]
+			Mjkplus3 = pipeparts.mkaudioamplify(pipeline, ones, factor)
+			factor = freqs[1] * (pow(freqs[0], 2) - pow(freqs[2 + j], 2)) * EPICS[1 + 2 * ((1 + num_stages) + 1 + k)] + freqs[0] * (pow(freqs[2 + j], 2) - pow(freqs[1], 2)) * EPICS[1 + 2 * (1 + k)] + freqs[2 + j] * (pow(freqs[1], 2) - pow(freqs[0], 2)) * EPICS[1 + 2 * ((2 + j) * (1 + num_stages) + 1 + k)]
+			Mjplus3k = pipeparts.mkaudioamplify(pipeline, ones, factor)
+			factor = 2.0 * numpy.pi * pow(freqs[1], 2) * (pow(freqs[0], 2) - pow(freqs[2 + j], 2)) * EPICS[2 * ((1 + num_stages) + 1 + k)] + 2.0 * numpy.pi * pow(freqs[0], 2) * (pow(freqs[2 + j], 2) - pow(freqs[1], 2)) * EPICS[2 * (1 + k)] + 2.0 * numpy.pi * pow(freqs[2 + j], 2) * (pow(freqs[1], 2) - pow(freqs[0], 2)) * EPICS[2 * ((2 + j) * (1 + num_stages) + 1 + k)]
+			Mjplus3kplus3 = pipeparts.mkaudioamplify(pipeline, ones, factor)
+			# Add these into the matrix
+			MV_matrix[(1 + j) * 2 * num_stages + k] = Mjk
+			MV_matrix[(1 + j) * 2 * num_stages + num_stages + k] = Mjkplus3
+			MV_matrix[(1 + num_stages + j) * 2 * num_stages + k] = Mjplus3k
+			MV_matrix[(1 + num_stages + j) * 2 * num_stages + num_stages + k] = Mjplus3kplus3
+
+	print("exact kappas 20")
+	# Now pass these to the matrix solver to find kappa_T, kappa_P, kappa_U, tau_T, tau_P, and tau_U.
+	MV_matrix = mkinterleave(pipeline, MV_matrix)
+	kappas = pipeparts.mkgeneric(pipeline, MV_matrix, "lal_matrixsolver")
+	kappas = list(mkdeinterleave(pipeline, kappas, 2 * num_stages))
+	for i in range(len(kappas)):
+		if i >= len(kappas) / 2:
+			kappas[i] = mkmultiplier(pipeline, [kappas[i], mkpow(pipeline, kappas[i - len(kappas) / 2], exponent = -1.0)])
+		kappas[i] = pipeparts.mktee(pipeline, kappas[i])
+
+	# Next, compute kappa_C.  This is going to take some work...
+	# Start by computing G_res at each frequency, defined in Eq. 5.2.30
+	print("exact kappas 25")
+	for n in range(2):
+		Gres_components = []
+		for j in range(num_stages):
+			kappajGresjatn = pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, kappas[j], matrix = [[EPICS[2 * (n * (1 + num_stages) + 1 + j)], EPICS[1 + 2 * (n * (1 + num_stages) + 1 + j)]]]))
+			i_omega_tau = pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, kappas[num_stages + j], matrix = [[0, 2.0 * numpy.pi * freqs[n]]]))
+			phase = pipeparts.mkgeneric(pipeline, i_omega_tau, "cexp")
+			Gres_components.append(mkmultiplier(pipeline, list_srcs(pipeline, kappajGresjatn, phase)))
+		Gres = pipeparts.mktee(pipeline, mkadder(pipeline, Gres_components))
+		Gresreal.append(pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkgeneric(pipeline, Gres, "creal"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0")))
+		Gresimag.append(pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkgeneric(pipeline, Gres, "cimag"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0")))
+
+	print("exact kappas 30")
+	# Next, let us find the H's, I's, Xi, and zeta, defined by Eqs. 5.2.48, 5.2.49, 5.2.50, and 5.2.58, respectively.
+	H1 = pipeparts.mktee(pipeline, pipeparts.mkaudioamplify(pipeline, mkadder(pipeline, list_srcs(pipeline, Gresreal[0], pipeparts.mkaudioamplify(pipeline, Yreal[0], -1.0))), pow(freqs[0], 2)))
+	print("exact kappas 31")
+	H2 = pipeparts.mktee(pipeline, pipeparts.mkaudioamplify(pipeline, mkadder(pipeline, list_srcs(pipeline, Gresreal[1], pipeparts.mkaudioamplify(pipeline, Yreal[1], -1.0))), pow(freqs[1], 2)))
+	print("exact kappas 32")
+	I1 = pipeparts.mktee(pipeline, pipeparts.mkaudioamplify(pipeline, mkadder(pipeline, list_srcs(pipeline, Gresimag[0], pipeparts.mkaudioamplify(pipeline, Yimag[0], -1.0))), pow(freqs[0], 2)))
+	print("exact kappas 33")
+	I2 = pipeparts.mktee(pipeline, pipeparts.mkaudioamplify(pipeline, mkadder(pipeline, list_srcs(pipeline, Gresimag[1], pipeparts.mkaudioamplify(pipeline, Yimag[1], -1.0))), pow(freqs[1], 2)))
+	print("exact kappas 34")
+	Xi = pipeparts.mktee(pipeline, mkadder(pipeline, list_srcs(pipeline, pipeparts.mkaudioamplify(pipeline, H1, pow(freqs[1], 2) / (pow(freqs[0], 2) - pow(freqs[1], 2))), pipeparts.mkaudioamplify(pipeline, H2, -pow(freqs[0], 2) / (pow(freqs[0], 2) - pow(freqs[1], 2))))))
+	print("exact kappas 35")
+	zeta = pipeparts.mktee(pipeline, mkadder(pipeline, list_srcs(pipeline, pipeparts.mkaudioamplify(pipeline, I1, 1.0 / freqs[0]), pipeparts.mkaudioamplify(pipeline, I2, -1.0 / freqs[1]))))
+
+	print("exact kappas 40")
+	# Now, we define the coefficients of the quartic equation a kC^4 + b kC^3 + c kC^2 + d kC + e = 0.
+	a = pipeparts.mkaudioamplify(pipeline, mkmultiplier(pipeline, [Xi, mkadder(pipeline, [H2, Xi]), mkpow(pipeline, zeta, exponent = 2.0)]), pow(freqs[1], 2))
+	print("exact kappas 41")
+	a = pipeparts.mktee(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, a, matrix = [[1.0, 0.0]])))
+	print("exact kappas 42")
+	b = mkadder(pipeline, [pipeparts.mkaudioamplify(pipeline, mkmultiplier(pipeline, [mkadder(pipeline, [H2, Xi]), zeta, I2]), freqs[1] * (pow(freqs[1], 2) - pow(freqs[0], 2))), pipeparts.mkaudioamplify(pipeline, mkmultiplier(pipeline, [mkadder(pipeline, [H2, pipeparts.mkaudioamplify(pipeline, Xi, 2.0)]), mkpow(pipeline, zeta, exponent = 2.0)]), pow(freqs[1], 4))])
+	print("exact kappas 43")
+	b = pipeparts.mktee(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, b, matrix = [[1.0, 0.0]])))
+	print("exact kappas 44")
+	c = mkadder(pipeline, [pipeparts.mkaudioamplify(pipeline, mkmultiplier(pipeline, [zeta, I2]), pow(freqs[1], 3) * (pow(freqs[1], 2) - pow(freqs[0], 2))), pipeparts.mkaudioamplify(pipeline, mkpow(pipeline, mkadder(pipeline, [H2, Xi]), exponent = 2.0), pow(pow(freqs[1], 2) - pow(freqs[0], 2), 2)), pipeparts.mkaudioamplify(pipeline, mkpow(pipeline, zeta, exponent = 2.0), pow(freqs[1], 6))])
+	print("exact kappas 45")
+	c = pipeparts.mktee(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, c, matrix = [[1.0, 0.0]])))
+	print("exact kappas 46")
+	d = pipeparts.mkaudioamplify(pipeline, mkadder(pipeline, [H2, Xi]), pow(freqs[1], 4) * pow(pow(freqs[1], 2) - pow(freqs[0], 2), 2))
+	d = pipeparts.mktee(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, d, matrix = [[1.0, 0.0]])))
+	e = pow(freqs[1], 4) * pow(pow(freqs[1], 2) - pow(freqs[0], 2), 2)
+
+	print("exact kappas 50")
+	# More parameters we need to define to solve the quartic...
+	Delta0 = pipeparts.mktee(pipeline, mkadder(pipeline, [mkpow(pipeline, c, exponent = 2.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [b, d]), -3.0, 0.0), complex_audioamplify(pipeline, a, 12 * e, 0.0)]))
+	print("exact kappas 51")
+	Delta1 = pipeparts.mktee(pipeline, mkadder(pipeline, [complex_audioamplify(pipeline, mkpow(pipeline, c, exponent = 3.0), 2.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [b, c, d]), -9.0, 0.0), complex_audioamplify(pipeline, mkpow(pipeline, b, exponent = 2.0), 27.0 * e, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [a, mkpow(pipeline, d, exponent = 2.0)]), 27.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [a, c]), -72.0 * e, 0.0)]))
+	print("exact kappas 52")
+	p = pipeparts.mktee(pipeline, mkmultiplier(pipeline, [mkadder(pipeline, [complex_audioamplify(pipeline, mkmultiplier(pipeline, [a, c]), 8.0, 0.0), complex_audioamplify(pipeline, mkpow(pipeline, b, exponent = 2.0), -3.0, 0.0)]), complex_audioamplify(pipeline, mkpow(pipeline, a, exponent = -2.0), 1.0 / 8.0, 0.0)]))
+	print("exact kappas 53")
+	q = mkmultiplier(pipeline, [mkadder(pipeline, [mkpow(pipeline, b, exponent = 3.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [a, b, c]), -4.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [mkpow(pipeline, a, exponent = 2.0), d]), 8.0, 0.0)]), complex_audioamplify(pipeline, mkpow(pipeline, a, exponent = -3.0), 1.0 / 8.0, 0.0)])
+	print("exact kappas 54")
+	Q0 = pipeparts.mktee(pipeline, mkpow(pipeline, complex_audioamplify(pipeline, mkadder(pipeline, [Delta1, mkpow(pipeline, mkadder(pipeline, [mkpow(pipeline, Delta1, exponent = 2.0), complex_audioamplify(pipeline, mkpow(pipeline, Delta0, exponent = 3.0), -4.0, 0.0)]), exponent = 0.5)]), 0.5, 0.0), exponent = 1.0 / 3.0))
+	print("exact kappas 55")
+	S0 = pipeparts.mktee(pipeline, mkpow(pipeline, mkadder(pipeline, [complex_audioamplify(pipeline, p, -1.0 / 6.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [mkpow(pipeline, a, exponent = -1.0), mkadder(pipeline, [Q0, complex_division(pipeline, Delta0, Q0)])]), 1.0 / 12.0, 0.0)]), exponent = 0.5))
+	print("exact kappas 56")
+	R0 = pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkgeneric(pipeline, mkadder(pipeline, [mkpow(pipeline, b, exponent = 3.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [d, mkpow(pipeline, a, exponent = 2.0)]), 8.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [a, b, c]), -4.0, 0.0)]), "creal"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0"))
+	print("exact kappas 57")
+	R0sign = pipeparts.mktee(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, mkmultiplier(pipeline, [R0, mkpow(pipeline, pipeparts.mkgeneric(pipeline, R0, "cabs"), exponent = -1.0)]), matrix = [[1.0, 0.0]])))
+
+	print("exact kappas 60")
+	# Finally, compute kappa_C
+	kappa_C = mkadder(pipeline, [complex_audioamplify(pipeline, complex_division(pipeline, b, a), -0.25, 0.0), mkmultiplier(pipeline, [R0sign, S0]), complex_audioamplify(pipeline, mkpow(pipeline, mkadder(pipeline, [complex_audioamplify(pipeline, mkpow(pipeline, S0, exponent = 2.0), -4.0, 0.0), complex_audioamplify(pipeline, p, -2.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [q, R0sign, mkpow(pipeline, S0, exponent = -1.0)]), -1.0, 0.0)]), exponent = 0.5), -0.5, 0.0)])
+	kappa_C = pipeparts.mktee(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkgeneric(pipeline, kappa_C, "creal"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0"))
+	kappas.append(kappa_C)
+
+	print("exact kappas 70")
+	# Next, compute f_cc using Eq. 5.2.45
+	f_cc = pipeparts.mktee(pipeline, pipeparts.mkaudioamplify(pipeline, mkpow(pipeline, mkmultiplier(pipeline, [kappa_C, mkadder(pipeline, [pipeparts.mkaudioamplify(pipeline, Gresimag[0], freqs[0]), pipeparts.mkaudioamplify(pipeline, Yimag[0], -freqs[0]), pipeparts.mkaudioamplify(pipeline, Gresimag[1], -freqs[1]), pipeparts.mkaudioamplify(pipeline, Yimag[1], freqs[1])])]), exponent = -1.0), pow(freqs[1], 2) - pow(freqs[0], 2)))
+	kappas.append(f_cc)
+
+	print("exact kappas 80")
+	# Compute f_s^2 using Eq. 5.2.41
+	f_s_squared = pipeparts.mktee(pipeline, pipeparts.mkaudioamplify(pipeline, mkmultiplier(pipeline, [kappa_C, mkadder(pipeline, [Yreal[1], pipeparts.mkaudioamplify(pipeline, Gresreal[1], -1.0), pipeparts.mkaudioamplify(pipeline, Yreal[0], -1.0), Gresreal[0]])]), 1.0 / (pow(freqs[1], -2) - pow(freqs[0], -2))))
+	kappas.append(f_s_squared)
+
+	print("exact kappas 90")
+	# Compute f_s / Q using Eq. 5.2.33
+	f_s_over_Q = pipeparts.mktee(pipeline, mkmultiplier(pipeline, [f_cc, pipeparts.mkgeneric(pipeline, mkadder(pipeline, [mkmultiplier(pipeline, [kappa_C, Yreal[0]]), pipeparts.mkaudioamplify(pipeline, f_s_squared, -pow(freqs[0], -2)), pipeparts.mkaudioamplify(pipeline, mkmultiplier(pipeline, [kappa_C, Gresreal[0]]), -1.0)]), "lal_add_constant", value = -1.0)]))
+	kappas.append(f_s_over_Q)
+
+	print("exact kappas done")
+	return kappas
 
 def compute_S_from_filters_file(pipeline, EP6R, EP6I, pcalfpcal2, derrfpcal2, EP7R, EP7I, ktst, EP8R, EP8I, kpu, EP9R, EP9I):
 
diff --git a/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py b/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py
index bf63ab5bfd644bdc518a71d1823cec69e9e79571..a0560e50b271983959180e0afa8fbcf62c866186 100644
--- a/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py
+++ b/gstlal-calibration/tests/check_calibration/pcal2darm_timeseries.py
@@ -80,6 +80,7 @@ parser.add_option("--phase-ranges", metavar = "list", type = str, default = "-6.
 parser.add_option("--labels", metavar = "list", type = str, help = "Comma-separated List of labels for each calibrated channel being tested. This is put in the plot legends and in the txt file names to distinguish them.")
 parser.add_option("--file-name-suffix", metavar = "name", type = str, default = "", help = "Suffix for naming unique file.")
 parser.add_option("--pcal-time-advance", metavar = "seconds", type = float, default = 0.0, help = "Time advance in seconds applied to the Pcal channel. Default = 0.0")
+parser.add_option("--show-stats", action = "store_true", help = "If set, plots will have averages (mean) and standard deviations shown in plot legends.")
 
 options, filenames = parser.parse_args()
 
@@ -315,7 +316,10 @@ for i in range(0, len(frequencies)):
 	if i == 0:
 		plt.figure(figsize = (25, 15))
 	plt.subplot(2, len(frequencies), i + 1)
-	plt.plot(times, magnitudes[0], colors[0], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu = %0.4f, \sigma = %0.4f]$' % (labels[0].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_'), numpy.mean(magnitudes[0]), numpy.std(magnitudes[0])))
+	if options.show_stats:
+		plt.plot(times, magnitudes[0], colors[0], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu = %0.4f, \sigma = %0.4f]$' % (labels[0].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_'), numpy.mean(magnitudes[0]), numpy.std(magnitudes[0])))
+	else:
+		plt.plot(times, magnitudes[0], colors[0], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s}$' % (labels[0].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_')))
 	plt.title(r'${\rm %s} \ \widetilde{\Delta L}_{\rm free} / \tilde{x}_{\rm pc} \  {\rm at \  %0.1f \  Hz}$' % ( ifo, frequencies[i]), fontsize = 32)
 	if i == 0:
 		plt.ylabel(r'${\rm Magnitude}$')
@@ -325,7 +329,10 @@ for i in range(0, len(frequencies)):
 	leg = plt.legend(fancybox = True, markerscale = 8.0 / markersize, numpoints = 3)
 	leg.get_frame().set_alpha(0.5)
 	plt.subplot(2, len(frequencies), len(frequencies) + i + 1)
-	plt.plot(times, phases[0], colors[0], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu = %0.3f^{\circ}, \sigma = %0.3f^{\circ}]$' % (labels[0].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_'), numpy.mean(phases[0]), numpy.std(phases[0])))
+	if options.show_stats:
+		plt.plot(times, phases[0], colors[0], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu = %0.3f^{\circ}, \sigma = %0.3f^{\circ}]$' % (labels[0].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_'), numpy.mean(phases[0]), numpy.std(phases[0])))
+	else:
+		plt.plot(times, phases[0], colors[0], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s}$' % (labels[0].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_')))
 	leg = plt.legend(fancybox = True, markerscale = 8.0 / markersize, numpoints = 3)
 	leg.get_frame().set_alpha(0.5)
 	if i == 0:
@@ -343,11 +350,17 @@ for i in range(0, len(frequencies)):
 			magnitudes[j].append(data[filter_time * k][1])
 			phases[j].append(data[filter_time * k][2])
 		plt.subplot(2, len(frequencies), i + 1)
-		plt.plot(times, magnitudes[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu = %0.4f, \sigma = %0.4f]$' % (labels[j].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_').replace(' ', '\ '), numpy.mean(magnitudes[j]), numpy.std(magnitudes[j])))
+		if options.show_stats:
+			plt.plot(times, magnitudes[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu = %0.4f, \sigma = %0.4f]$' % (labels[j].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_').replace(' ', '\ '), numpy.mean(magnitudes[j]), numpy.std(magnitudes[j])))
+		else:
+			plt.plot(times, magnitudes[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s}$' % (labels[j].replace(':', '{:}').replace('-', '\mbox{-}').replace('_', '\_').replace(' ', '\ ')))
 		leg = plt.legend(fancybox = True, markerscale = 8.0 / markersize, numpoints = 3)
 		leg.get_frame().set_alpha(0.5)
 		plt.subplot(2, len(frequencies), len(frequencies) + i + 1)
-		plt.plot(times, phases[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu = %0.3f^{\circ}, \sigma = %0.3f^{\circ}]$' % (labels[j].replace(':', '{:}').replace('-' , '\mbox{-}').replace('_', '\_').replace(' ', '\ '), numpy.mean(phases[j]), numpy.std(phases[j])))
+		if options.show_stats:
+			plt.plot(times, phases[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s} \ [\mu = %0.3f^{\circ}, \sigma = %0.3f^{\circ}]$' % (labels[j].replace(':', '{:}').replace('-' , '\mbox{-}').replace('_', '\_').replace(' ', '\ '), numpy.mean(phases[j]), numpy.std(phases[j])))
+		else:
+			plt.plot(times, phases[j], colors[j % 6], linestyle = 'None', marker = '.', markersize = markersize, label = r'${\rm %s}$' % (labels[j].replace(':', '{:}').replace('-' , '\mbox{-}').replace('_', '\_').replace(' ', '\ ')))
 		leg = plt.legend(fancybox = True, markerscale = 8.0 / markersize, numpoints = 3)
 		leg.get_frame().set_alpha(0.5)
 plt.savefig("%s_deltal_over_pcal%s_%d-%d.png" % (ifo, options.file_name_suffix, int(t_start), int(dur_in_seconds)))
diff --git a/gstlal-calibration/tests/test_common.py b/gstlal-calibration/tests/test_common.py
index b0f2c7bc2bd7269f59c911a14f890471b709a462..48cec7604a32a8f8ddc931fd21ded6d0e2891a48 100644
--- a/gstlal-calibration/tests/test_common.py
+++ b/gstlal-calibration/tests/test_common.py
@@ -58,12 +58,12 @@ else:
 #
 
 
-def complex_test_src(pipeline, buffer_length = 1.0, rate = 2048, width = 64, test_duration = 10.0, wave = 5, freq = 0, is_live = False, verbose = True):
+def complex_test_src(pipeline, buffer_length = 1.0, rate = 2048, width = 64, test_duration = 10.0, wave = 5, freq = 0, is_live = False, verbose = True, src_suffix = ""):
 	head = pipeparts.mkaudiotestsrc(pipeline, wave = wave, freq = freq, samplesperbuffer = int(buffer_length * rate), volume = 1, num_buffers = int(test_duration / buffer_length), is_live = is_live)
-	head = pipeparts.mkcapsfilter(pipeline, head, "audio/x-raw, format=Z%d%s, rate=%d, channels=2" % (width, BYTE_ORDER, rate))
+	head = pipeparts.mkcapsfilter(pipeline, head, "audio/x-raw, format=F%d%s, rate=%d, channels=2" % (width, BYTE_ORDER, rate))
 	head = pipeparts.mktogglecomplex(pipeline, head)
 	if verbose:
-		head = pipeparts.mkprogressreport(pipeline, head, "src")
+		head = pipeparts.mkprogressreport(pipeline, head, "src%s" % src_suffix)
 	return head
 
 def int_test_src(pipeline, buffer_length = 1.0, rate = 2048, width = 64, channels = 1, test_duration = 10.0, wave = 5, freq = 0, is_live = False, verbose = True):