Commit ed911fcc authored by Aaron Viets's avatar Aaron Viets
Browse files

gstlal_compute_strain: Exact kappas solution runs now... but results not sensical yet

parent 6ee2014b
Pipeline #74569 failed with stages
in 1 minute and 14 seconds
...@@ -1425,23 +1425,23 @@ if compute_exact_kappas and (compute_kappatst or compute_kappapum or compute_kap ...@@ -1425,23 +1425,23 @@ if compute_exact_kappas and (compute_kappatst or compute_kappapum or compute_kap
pipeparts.mkfakesink(pipeline, tau_uim) pipeparts.mkfakesink(pipeline, tau_uim)
if compute_kappac: if compute_kappac:
smooth_kctee = smooth_kc_nogate = kc smooth_kctee = smooth_kc_nogate = pipeparts.mktee(pipeline, kc)
else: else:
pipeparts.mkfakesink(pipeline, kc) pipeparts.mkfakesink(pipeline, kc)
if compute_fcc: if compute_fcc:
smooth_fcctee = smooth_fcc_nogate = fcc smooth_fcctee = smooth_fcc_nogate = pipeparts.mktee(pipeline, fcc)
else: else:
pipeparts.mkfakesink(pipeline, fcc) pipeparts.mkfakesink(pipeline, fcc)
if compute_fs: if compute_fs:
smooth_fs_squared_nogate = smooth_fs_squared = fs_squared smooth_fs_squared_nogate = smooth_fs_squared = pipeparts.mktee(pipeline, fs_squared)
else: else:
pipeparts.mkfakesink(pipeline, fs_squared) pipeparts.mkfakesink(pipeline, fs_squared)
if compute_srcq: 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. # 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)])) smooth_srcQ_inv = 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, smooth_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 # 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 # 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 # antispring, the computed Q^(-1) is real. To get a sensible real value either
......
...@@ -459,11 +459,11 @@ static gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, GstCaps *outc ...@@ -459,11 +459,11 @@ static gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, GstCaps *outc
if(!strcmp(name, GST_AUDIO_NE(F32))) { if(!strcmp(name, GST_AUDIO_NE(F32))) {
element->data_type = GSTLAL_SENSINGTDCFS_FLOAT; element->data_type = GSTLAL_SENSINGTDCFS_FLOAT;
g_assert_cmpuint(unit_size_out, ==, 4 * (guint) channels_out); g_assert_cmpuint(unit_size_out, ==, 4 * (guint) channels_out);
g_assert_cmpuint(unit_size_in, ==, 4 * (guint) channels_in); g_assert_cmpuint(unit_size_in, ==, 8 * (guint) channels_in);
} else if(!strcmp(name, GST_AUDIO_NE(F64))) { } else if(!strcmp(name, GST_AUDIO_NE(F64))) {
element->data_type = GSTLAL_SENSINGTDCFS_DOUBLE; element->data_type = GSTLAL_SENSINGTDCFS_DOUBLE;
g_assert_cmpuint(unit_size_out, ==, 8 * (guint) channels_out); g_assert_cmpuint(unit_size_out, ==, 8 * (guint) channels_out);
g_assert_cmpuint(unit_size_in, ==, 8 * (guint) channels_in); g_assert_cmpuint(unit_size_in, ==,16 * (guint) channels_in);
} else } else
g_assert_not_reached(); g_assert_not_reached();
......
...@@ -975,26 +975,25 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate): ...@@ -975,26 +975,25 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate):
# ratios X[i] = injection(f_i) / d_err(f_i) for each calibration line frequency. # ratios X[i] = injection(f_i) / d_err(f_i) for each calibration line frequency.
# #
print("starting exact kappas")
kappas = [] kappas = []
num_lines = len(freqs) num_lines = len(freqs)
num_stages = num_lines - 2 # Stages of actuation (currently 3) num_stages = num_lines - 2 # Stages of actuation (currently 3)
MV_matrix = list(numpy.zeros(2 * num_stages * (2 * num_stages + 1))) MV_matrix = list(numpy.zeros(2 * num_stages * (2 * num_stages + 1)))
Y = []
Yreal = [] Yreal = []
Yimag = [] Yimag = []
CAX = [] CAX = []
CAXreal = [] CAXreal = []
CAXimag = [] CAXimag = []
Gresreal = [] Gres = []
Gresimag = []
kappas = [] kappas = []
for i in range(num_lines): for i in range(num_lines):
if i < 2: if i < 2:
# Then it's a Pcal line # 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])) Y.append(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, mkcapsfiltersetter(pipeline, pipeparts.mkgeneric(pipeline, Y, "creal"), "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate, name = "capsfilter_Yreal_%d" % i))) Yreal.append(pipeparts.mktee(pipeline, mkcapsfiltersetter(pipeline, pipeparts.mkgeneric(pipeline, Y[i], "creal"), "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate, name = "capsfilter_Yreal_%d" % i)))
Yimag.append(pipeparts.mktee(pipeline, mkcapsfiltersetter(pipeline, pipeparts.mkgeneric(pipeline, Y, "cimag"), "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate, name = "capsfilter_Yimag_%d" % i))) Yimag.append(pipeparts.mktee(pipeline, mkcapsfiltersetter(pipeline, pipeparts.mkgeneric(pipeline, Y[i], "cimag"), "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate, name = "capsfilter_Yimag_%d" % i)))
else: else:
# It's an actuator line # 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]))) CAX.append(pipeparts.mktee(pipeline, complex_audioamplify(pipeline, X[i], EPICS[2 * (1 + num_stages) * i], EPICS[2 * (1 + num_stages) * i + 1])))
...@@ -1018,7 +1017,6 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate): ...@@ -1018,7 +1017,6 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate):
# Many of the elements are constant, so make a stream of ones to multiply # Many of the elements are constant, so make a stream of ones to multiply
if num_stages > 1: if num_stages > 1:
ones = pipeparts.mktee(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, Yreal[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 0.0)) ones = pipeparts.mktee(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, Yreal[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 0.0))
print("exact kappas 10")
for j in range(num_stages): for j in range(num_stages):
# Time-dependent matrix elements # Time-dependent matrix elements
factor = pow(freqs[0], -2) - pow(freqs[1], -2) factor = pow(freqs[0], -2) - pow(freqs[1], -2)
...@@ -1057,7 +1055,6 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate): ...@@ -1057,7 +1055,6 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate):
MV_matrix[(1 + num_stages + j) * 2 * num_stages + k] = Mjplus3k MV_matrix[(1 + num_stages + j) * 2 * num_stages + k] = Mjplus3k
MV_matrix[(1 + num_stages + j) * 2 * num_stages + num_stages + k] = Mjplus3kplus3 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. # 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) MV_matrix = mkinterleave(pipeline, MV_matrix)
MV_matrix = pipeparts.mkcapsfilter(pipeline, MV_matrix, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=%d" % (rate, (2 * num_stages) * (2 * num_stages + 1))) MV_matrix = pipeparts.mkcapsfilter(pipeline, MV_matrix, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=%d" % (rate, (2 * num_stages) * (2 * num_stages + 1)))
...@@ -1071,7 +1068,6 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate): ...@@ -1071,7 +1068,6 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate):
# Next, compute kappa_C. This is going to take some work... # 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 # Start by computing G_res at each frequency, defined in Eq. 5.2.30
print("exact kappas 25")
for n in range(2): for n in range(2):
Gres_components = [] Gres_components = []
for j in range(num_stages): for j in range(num_stages):
...@@ -1080,82 +1076,14 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate): ...@@ -1080,82 +1076,14 @@ def compute_exact_kappas_from_filters_file(pipeline, X, freqs, EPICS, rate):
i_omega_tau = mkcapsfiltersetter(pipeline, i_omega_tau, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate) i_omega_tau = mkcapsfiltersetter(pipeline, i_omega_tau, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)
phase = pipeparts.mkgeneric(pipeline, i_omega_tau, "cexp") phase = pipeparts.mkgeneric(pipeline, i_omega_tau, "cexp")
Gres_components.append(mkmultiplier(pipeline, list_srcs(pipeline, kappajGresjatn, phase))) Gres_components.append(mkmultiplier(pipeline, list_srcs(pipeline, kappajGresjatn, phase)))
Gres = pipeparts.mktee(pipeline, mkadder(pipeline, Gres_components)) Gres.append(mkadder(pipeline, Gres_components))
Gresreal.append(pipeparts.mktee(pipeline, mkcapsfiltersetter(pipeline, pipeparts.mkgeneric(pipeline, Gres, "creal"), "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)))
Gresimag.append(pipeparts.mktee(pipeline, mkcapsfiltersetter(pipeline, pipeparts.mkgeneric(pipeline, Gres, "cimag"), "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate))) sensing_inputs = mkinterleave(pipeline, Gres + Y, complex_data = True)
sensing_outputs = pipeparts.mkgeneric(pipeline, sensing_inputs, "lal_sensingtdcfs", sensing_model = 0, freq1 = freqs[0], freq2 = freqs[1])
print("exact kappas 30") sensing_outputs = list(mkdeinterleave(pipeline, sensing_outputs, 4))
# 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.mkaudioamplify(pipeline, mkadder(pipeline, list_srcs(pipeline, pipeparts.mkcapsfilter(pipeline, Gresreal[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkcapsfilter(pipeline, Yreal[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -1.0))), pow(freqs[0], 2)) kappas += sensing_outputs
print("exact kappas 31")
H2 = pipeparts.mktee(pipeline, pipeparts.mkaudioamplify(pipeline, mkadder(pipeline, list_srcs(pipeline, pipeparts.mkcapsfilter(pipeline, Gresreal[1], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, pipeparts.mkcapsfilter(pipeline, Yreal[1], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -1.0))), pow(freqs[1], 2)))
print("exact kappas 32")
I1 = pipeparts.mkaudioamplify(pipeline, mkadder(pipeline, list_srcs(pipeline, pipeparts.mkcapsfilter(pipeline, Gresimag[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, Yimag[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -1.0))), pow(freqs[0], 2))
print("exact kappas 33")
I2 = pipeparts.mktee(pipeline, pipeparts.mkaudioamplify(pipeline, mkadder(pipeline, list_srcs(pipeline, pipeparts.mkcapsfilter(pipeline, Gresimag[1], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, Yimag[1], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -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, pipeparts.mkcapsfilter(pipeline, H2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -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, pipeparts.mkcapsfilter(pipeline, I2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -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, [pipeparts.mkcapsfilter(pipeline, Xi, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), mkadder(pipeline, [pipeparts.mkcapsfilter(pipeline, H2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, Xi, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, zeta, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), 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, [pipeparts.mkcapsfilter(pipeline, H2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, Xi, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), pipeparts.mkcapsfilter(pipeline, zeta, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, I2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), freqs[1] * (pow(freqs[1], 2) - pow(freqs[0], 2))), pipeparts.mkaudioamplify(pipeline, mkmultiplier(pipeline, [mkadder(pipeline, [pipeparts.mkcapsfilter(pipeline, H2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, Xi, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), 2.0)]), mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, zeta, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), 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, [pipeparts.mkcapsfilter(pipeline, zeta, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, I2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), pow(freqs[1], 3) * (pow(freqs[1], 2) - pow(freqs[0], 2))), pipeparts.mkaudioamplify(pipeline, mkpow(pipeline, mkadder(pipeline, [pipeparts.mkcapsfilter(pipeline, H2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, Xi, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), exponent = 2.0), pow(pow(freqs[1], 2) - pow(freqs[0], 2), 2)), pipeparts.mkaudioamplify(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, zeta, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), 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, [pipeparts.mkcapsfilter(pipeline, H2, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, Xi, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), 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, pipeparts.mkcapsfilter(pipeline, c, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 2.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, d, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), -3.0, 0.0), complex_audioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), 12 * e, 0.0)]))
print("exact kappas 51")
Delta1 = pipeparts.mktee(pipeline, mkadder(pipeline, [complex_audioamplify(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, c, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 3.0), 2.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, c, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, d, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), -9.0, 0.0), complex_audioamplify(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 2.0), 27.0 * e, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, d, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 2.0)]), 27.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, c, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), -72.0 * e, 0.0)]))
print("exact kappas 52")
p = pipeparts.mktee(pipeline, mkmultiplier(pipeline, [mkadder(pipeline, [complex_audioamplify(pipeline, mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, c, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), 8.0, 0.0), complex_audioamplify(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 2.0), -3.0, 0.0)]), complex_audioamplify(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = -2.0), 1.0 / 8.0, 0.0)]))
print("exact kappas 53")
q = mkmultiplier(pipeline, [mkadder(pipeline, [mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 3.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, c, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), -4.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 2.0), pipeparts.mkcapsfilter(pipeline, d, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), 8.0, 0.0)]), complex_audioamplify(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = -3.0), 1.0 / 8.0, 0.0)])
print("exact kappas 54")
Q0 = pipeparts.mktee(pipeline, mkpow(pipeline, complex_audioamplify(pipeline, mkadder(pipeline, [pipeparts.mkcapsfilter(pipeline, Delta1, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), mkpow(pipeline, mkadder(pipeline, [mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, Delta1, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 2.0), complex_audioamplify(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, Delta0, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), 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, pipeparts.mkcapsfilter(pipeline, p, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -1.0 / 6.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = -1.0), mkadder(pipeline, [pipeparts.mkcapsfilter(pipeline, Q0, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), complex_division(pipeline, pipeparts.mkcapsfilter(pipeline, Delta0, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, Q0, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate))])]), 1.0 / 12.0, 0.0)]), exponent = 0.5))
print("exact kappas 56")
R0 = pipeparts.mktee(pipeline, mkcapsfiltersetter(pipeline, pipeparts.mkgeneric(pipeline, mkadder(pipeline, [mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 3.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, d, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 2.0)]), 8.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, c, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), -4.0, 0.0)]), "creal"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0,channels=1"))
print("exact kappas 57")
R0sign = pipeparts.mktee(pipeline, pipeparts.mktogglecomplex(pipeline, pipeparts.mkmatrixmixer(pipeline, mkmultiplier(pipeline, [pipeparts.mkgeneric(pipeline, pipeparts.mkcapsfilter(pipeline, R0, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), "creal"), mkpow(pipeline, pipeparts.mkgeneric(pipeline, pipeparts.mkcapsfilter(pipeline, R0, "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), "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, pipeparts.mkcapsfilter(pipeline, b, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, a, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)), -0.25, 0.0), mkmultiplier(pipeline, [pipeparts.mkcapsfilter(pipeline, R0sign, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkcapsfilter(pipeline, S0, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), complex_audioamplify(pipeline, mkpow(pipeline, mkadder(pipeline, [complex_audioamplify(pipeline, mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, S0, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = 2.0), -4.0, 0.0), complex_audioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, p, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -2.0, 0.0), complex_audioamplify(pipeline, mkmultiplier(pipeline, [q, pipeparts.mkcapsfilter(pipeline, R0sign, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), mkpow(pipeline, pipeparts.mkcapsfilter(pipeline, S0, "audio/x-raw,format=Z128LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), exponent = -1.0)]), -1.0, 0.0)]), exponent = 0.5), -0.5, 0.0)])
kappa_C = pipeparts.mktee(pipeline, mkcapsfiltersetter(pipeline, pipeparts.mkgeneric(pipeline, kappa_C, "creal"), "audio/x-raw,format=F64LE,channel-mask=(bitmask)0x0,channels=1"))
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, pipeparts.mkcapsfilter(pipeline, Gresimag[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), freqs[0]), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, Yimag[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -freqs[0]), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, Gresimag[1], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -freqs[1]), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, Yimag[1], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), 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, [pipeparts.mkcapsfilter(pipeline, Yreal[1], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, Gresreal[1], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -1.0), pipeparts.mkaudioamplify(pipeline, pipeparts.mkcapsfilter(pipeline, Yreal[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate), -1.0), pipeparts.mkcapsfilter(pipeline, Gresreal[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)])]), 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, pipeparts.mkcapsfilter(pipeline, Yreal[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), pipeparts.mkaudioamplify(pipeline, f_s_squared, -pow(freqs[0], -2)), pipeparts.mkaudioamplify(pipeline, mkmultiplier(pipeline, [kappa_C, pipeparts.mkcapsfilter(pipeline, Gresreal[0], "audio/x-raw,format=F64LE,rate=%d,channel-mask=(bitmask)0x0,channels=1" % rate)]), -1.0)]), "lal_add_constant", value = -1.0)]))
kappas.append(f_s_over_Q)
print("exact kappas done")
return kappas return kappas
def compute_S_from_filters_file(pipeline, EP6R, EP6I, pcalfpcal2, derrfpcal2, EP7R, EP7I, ktst, EP8R, EP8I, kpu, EP9R, EP9I): def compute_S_from_filters_file(pipeline, EP6R, EP6I, pcalfpcal2, derrfpcal2, EP7R, EP7I, ktst, EP8R, EP8I, kpu, EP9R, EP9I):
......
...@@ -8,9 +8,9 @@ IFO = H ...@@ -8,9 +8,9 @@ IFO = H
# determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14) # determines where to look for filters files (e.g., O1, O2, O3, ER10, ER13, ER14, PreER10, PreER13, PreER14)
OBSRUN = O3 OBSRUN = O3
START = $(shell echo 1246226240 - 256 | bc) START = $(shell echo 1246224240 - 256 | bc)
#1238997618 #1238997618
END = $(shell echo 1246226240 + 120 + 256 | bc) END = $(shell echo 1246225240 + 256 | bc)
#1240876818 #1240876818
SHMRUNTIME = 600 SHMRUNTIME = 600
# How much time does the calibration need to settle at the start and end? # How much time does the calibration need to settle at the start and end?
...@@ -21,6 +21,7 @@ GDSCONFIGS = Filters/O3/GDSFilters/H1GDS_1239476409_testAllCorrections.ini ...@@ -21,6 +21,7 @@ GDSCONFIGS = Filters/O3/GDSFilters/H1GDS_1239476409_testAllCorrections.ini
LOWLATENCYCONFIGS = Filters/O3/GDSFilters/H1GDS_lowLatency_1239476409.ini LOWLATENCYCONFIGS = Filters/O3/GDSFilters/H1GDS_lowLatency_1239476409.ini
LOWLATENCYASYMMETRICCONFIGS = Filters/O3/GDSFilters/H1GDS_lowLatency_asymmetric_1239476409.ini LOWLATENCYASYMMETRICCONFIGS = Filters/O3/GDSFilters/H1GDS_lowLatency_asymmetric_1239476409.ini
DCSCONFIGS = Filters/O3/GDSFilters/H1DCS_C01_1239472998_test.ini DCSCONFIGS = Filters/O3/GDSFilters/H1DCS_C01_1239472998_test.ini
EXACTKAPPASCONFIGS = Filters/O3/GDSFilters/H1DCS_exactKappas_1239472998.ini
DCSLINESCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_CleaningLines.ini DCSLINESCONFIGS = ../../config_files/O2/H1/tests/H1DCS_AllCorrections_CleaningLines.ini
#../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini #../../config_files/O2/H1/tests/H1DCS_AllCorrections_Cleaning.ini
DCSFCCCONFIGS = ../../config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini DCSFCCCONFIGS = ../../config_files/O2/H1/tests/H1DCS_FreqIndepAndFccCorrections_Cleaning.ini
...@@ -31,7 +32,7 @@ GDSOLDCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_old.ini ...@@ -31,7 +32,7 @@ GDSOLDCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_old.ini
GDSBETTERCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_better.ini GDSBETTERCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_better.ini
GDSBESTCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_best.ini GDSBESTCONFIGS = Filters/ER14/GDSFilters/L1GDS_1235491416_best.ini
all: pcal_DCS_transfer_functions all: $(IFO)1_hoft_DCS_EXACTKAPPAS_frames.cache
############################################### ###############################################
### These commands should change less often ### ### These commands should change less often ###
...@@ -87,6 +88,10 @@ $(IFO)1_hoft_DCS_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir ...@@ -87,6 +88,10 @@ $(IFO)1_hoft_DCS_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir
GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/DCS/ --frame-duration=64 --frames-per-file=1 --wings=256 --config-file $(DCSCONFIGS) GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/DCS/ --frame-duration=64 --frames-per-file=1 --wings=256 --config-file $(DCSCONFIGS)
ls Frames/$(OBSRUN)/$(IFO)1/DCS/$(IFO)-$(IFO)1DCS-*.gwf | lalapps_path2cache > $@ ls Frames/$(OBSRUN)/$(IFO)1/DCS/$(IFO)-$(IFO)1DCS-*.gwf | lalapps_path2cache > $@
$(IFO)1_hoft_DCS_EXACTKAPPAS_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir
GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/DCS/ --frame-duration=64 --frames-per-file=1 --wings=256 --config-file $(EXACTKAPPASCONFIGS)
ls Frames/$(OBSRUN)/$(IFO)1/DCS/$(IFO)-$(IFO)1DCS_exactKappas-*.gwf | lalapps_path2cache > $@
# In case we want to compare one calibration to another... # In case we want to compare one calibration to another...
$(IFO)1_hoft_GDS_TEST_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir $(IFO)1_hoft_GDS_TEST_frames.cache: $(IFO)1_easy_raw_frames.cache filters framesdir
GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/GDS/ --frame-duration=64 --frames-per-file=1 --wings=0 --config-file $(GDSTESTCONFIGS) GST_DEBUG=3 gstlal_compute_strain --gps-start-time $(START) --gps-end-time $(END) --frame-cache $(IFO)1_easy_raw_frames.cache --output-path Frames/$(OBSRUN)/$(IFO)1/GDS/ --frame-duration=64 --frames-per-file=1 --wings=0 --config-file $(GDSTESTCONFIGS)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment