diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py index 8328c71895812d5dbba08300db6cef0cbeb3ba17..fa420b48b3664c19b48adb6b4885a7b00360899f 100755 --- a/lisainstrument/instrument.py +++ b/lisainstrument/instrument.py @@ -87,8 +87,8 @@ class Instrument: # Optical pathlength noises backlink_asds=3E-12, backlink_fknees=2E-3, testmass_asds=2.4E-15, testmass_fknees=0.4E-3, oms_asds=(6.35E-12, 1.25E-11, 1.42E-12, 3.38E-12, 3.32E-12, 7.90E-12), oms_fknees=2E-3, - # TCB synchronization - sync_asds=0.42, + # MOC time correlation + moc_time_correlation_asds=0.42, # Tilt-to-length (TTL) ttl_coeffs='default', sc_jitter_asds=(5E-9, 5E-9, 5E-9), sc_jitter_fknees=(8E-4, 8E-4, 8E-4), @@ -147,16 +147,17 @@ class Instrument: clock_freqoffsets: dictionary of clock frequency offsets [s^-1], or 'default' clock_freqlindrifts: dictionary of clock frequency linear drifts [s^-2], or 'default' clock_freqquaddrifts: dictionary of clock frequency quadratic drifts [s^-2], or 'default' - clockinv_tolerance: convergence tolerance for timer deviation inversion [s] - clockinv_maxiter: maximum number of iterations for timer deviation inversion + clockinv_tolerance: convergence tolerance for clock noise inversion [s] + clockinv_maxiter: maximum number of iterations for clock noise inversion backlink_asds: dictionary of amplitude spectral densities for backlink noise [m/sqrt(Hz)] backlink_fknees: dictionary of cutoff frequencied for backlink noise [Hz] testmass_asds: dictionary of amplitude spectral densities for test-mass noise [ms^(-2)/sqrt(Hz)] testmass_fknees: dictionary of cutoff frequencies for test-mass noise [Hz] oms_asds: tuple of dictionaries of amplitude spectral densities for OMS noise [m/sqrt(Hz)], ordered as (isi_carrier, isi_usb, tmi_carrier, tmi_usb, rfi_carrier, rfi_usb) - sync_asds: dictionary of amplitude spectral densities for TCB synchronization noise [s/sqrt(Hz)], - the default ASD seems rather high, this is due to the small sampling rate (default 1 / 86400s) + moc_time_correlation_asds: dictionary of amplitude spectral densities for MOC time + correlation noise [s/sqrt(Hz)], the default ASD seems rather high, this is due + to the small sampling rate (default 1 / 86400s) oms_fknees: dictionary of cutoff frequencies for OMS noise ttl_coeffs: tuple (local_phi, distant_phi, local_eta, distant_eta) of dictionaries of tilt-to-length coefficients on each MOSA [m/rad], 'default' for a default set of @@ -318,8 +319,8 @@ class Instrument: else: self.clock_freqquaddrifts = ForEachSC(clock_freqquaddrifts) - # TCB synchronization - self.sync_asds = ForEachSC(sync_asds) + # MOC time correlation + self.moc_time_correlation_asds = ForEachSC(moc_time_correlation_asds) # Clock-noise inversion self.clockinv_tolerance = float(clockinv_tolerance) @@ -588,7 +589,7 @@ class Instrument: '13': 8.33159404, '32': 8.30446786, '21': 8.33159402, }) self.d_pprs = ForEachMOSA(0) - self.tps_proper_time_deviations = ForEachSC(0) + self.tps_wrt_tcb = ForEachSC(0) self.orbit_dataset = None elif isinstance(orbits, str): logger.info("Using orbit file '%s'", orbits) @@ -615,7 +616,7 @@ class Instrument: self.d_pprs = self.pprs.transformed(lambda _, x: 0 if np.isscalar(x) else np.gradient(x, self.physics_dt) ) - self.tps_proper_time_deviations = ForEachSC(0) + self.tps_wrt_tcb = ForEachSC(0) self.orbit_dataset = None def init_orbits_file_1_0(self, orbitf): @@ -643,7 +644,7 @@ class Instrument: raise ValueError(f"invalid orbit dataset '{self.orbit_dataset}'") return InterpolatedUnivariateSpline(times, values, k=5, ext='raise') - def tps_proper_time_deviations(sc): + def tps_wrt_tcb(sc): if self.orbit_dataset == 'tcb/ltt': return lambda x: 0 if self.orbit_dataset == 'tps/ppr': @@ -664,8 +665,8 @@ class Instrument: self.pprs = ForEachMOSA(lambda mosa: pprs(mosa)(self.physics_t)) logger.debug("Interpolating proper pseudo-range derivatives") self.d_pprs = ForEachMOSA(lambda mosa: d_pprs(mosa)(self.physics_t)) - logger.debug("Interpolating proper time deviation from TCB") - self.tps_proper_time_deviations = ForEachSC(lambda sc: tps_proper_time_deviations(sc)(self.telemetry_t)) + logger.debug("Interpolating TPSs with respect to TCB") + self.tps_wrt_tcb = ForEachSC(lambda sc: tps_wrt_tcb(sc)(self.physics_t_withinitial)) except ValueError as error: logger.error("Missing orbit information at \n%s", self.physics_t) raise ValueError("missing orbit information, use longer orbit file or adjust sampling") from error @@ -690,12 +691,12 @@ class Instrument: logger.debug("Interpolating proper pseudo-range derivatives") dataset = orbitf['tcb/d_ltt'] if self.orbit_dataset == 'tcb/ltt' else orbitf['tps/d_ppr'] self.d_pprs = ForEachMOSA(lambda mosa: interpolate(dataset[:, link_index[mosa]], self.physics_t)) - logger.debug("Interpolating proper time deviation from TCB") + logger.debug("Interpolating TPSs with respect to TCB") if self.orbit_dataset == 'tcb/ltt': - self.tps_proper_time_deviations = ForEachSC(lambda sc: 0) + self.tps_wrt_tcb = ForEachSC(lambda sc: 0) else: dataset = orbitf['tcb/delta_tau'] - self.tps_proper_time_deviations = ForEachSC(lambda sc: interpolate(dataset[:, sc_index[sc]], self.physics_t)) + self.tps_wrt_tcb = ForEachSC(lambda sc: interpolate(dataset[:, sc_index[sc]], self.physics_t)) except ValueError as error: logger.error("Missing orbit information at \n%s", self.physics_t) raise ValueError("missing orbit information, use longer orbit file or adjust sampling") from error @@ -850,10 +851,7 @@ class Instrument: self.sync_asds = ForEachSC(0) def disable_clock_noises(self): - """Turn off all imperfections on clocks. - - This includes offsets, and frequency offsets and deviations. - """ + """Turn off all imperfections on clocks.""" self.clock_asds = ForEachSC(0) self.clock_offsets = ForEachSC(0) self.clock_freqoffsets = ForEachSC(0) @@ -902,37 +900,41 @@ class Instrument: self.simulate_noises() - logger.debug("Computing local timer deviations") - t = self.physics_et_withinitial - self.local_timer_deviations_withinitial = \ + ## MOC time correlations + + logger.debug("Computing local THE with respect to TPS") + t = self.physics_et + self.the_wrt_tps_local = \ self.clock_offsets \ + self.clock_freqoffsets * t \ + self.clock_freqlindrifts * t**2 / 2 \ + self.clock_freqquaddrifts * t**3 / 3 \ - + ForEachSC(lambda sc: - cumulative_trapezoid(np.broadcast_to( - self.clock_noise_fluctuations_withinitial[sc], - self.physics_size + self.initial_telemetry_physics_size), - dx=self.physics_dt, initial=0) - ) - - self.local_timer_deviations = self.local_timer_deviations_withinitial.transformed( - lambda _, x: x if np.isscalar(x) else x[self.initial_telemetry_physics_size:] - ) - - ## Timer deviations from TCB + + self.integrated_clock_noise_fluctuations + logger.debug("Computing THE with respect to TCB") + t = self.physics_et_withinitial + self.the_wrt_tcb_withinitial = \ + self.tps_wrt_tcb \ + + self.clock_offsets \ + + self.clock_freqoffsets * (t + self.tps_wrt_tcb) \ + + self.clock_freqlindrifts * (t + self.tps_wrt_tcb)**2 / 2 \ + + self.clock_freqquaddrifts * (t + self.tps_wrt_tcb)**3 / 3 \ + + self.tps_wrt_tcb * self.clock_noise_fluctuations_withinitial \ + + self.integrated_clock_noise_fluctuations_withinitial + + logger.debug("Computing MOC time correlations") physics_to_telemetry = lambda _, x: x[::self.physics_upsampling * self.telemetry_downsampling] - self.tcb_timer_deviations = self.tcb_sync_noises \ - + self.local_timer_deviations_withinitial.transformed(physics_to_telemetry) + self.moc_time_correlations = self.moc_time_correlation_noises \ + + self.the_wrt_tcb_withinitial.transformed(physics_to_telemetry) ## TDIR tone + logger.debug("Forming TDIR assistance tones") self.tdir_tones = ForEachMOSA(lambda mosa: 0 if self.tdir_tone_amplitudes[mosa] == 0 \ else self.tdir_tone_amplitudes[mosa] * np.sin( 2 * np.pi * self.tdir_tone_frequencies[mosa] - * (self.physics_et + self.local_timer_deviations[mosa[0]]) + * (self.physics_et + self.the_wrt_tps_local[mosa[0]]) + self.tdir_tone_initial_phases[mosa] ) ) @@ -993,9 +995,9 @@ class Instrument: - (self.central_freq + self.local_usb_offsets) * self.gws \ - (self.central_freq + self.local_usb_offsets) * self.local_ttls / c - logger.debug("Propagating timer deviations to distant MOSAs") - self.distant_timer_deviations = \ - self.local_timer_deviations.for_each_mosa().distant() \ + logger.debug("Propagating local THEs with respect to TPS to distant MOSAs") + self.the_wrt_tps_distant = \ + self.the_wrt_tps_local.for_each_mosa().distant() \ .transformed(lambda mosa, x: self.interpolate(x, -self.pprs[mosa]) - self.pprs[mosa] ) @@ -1083,9 +1085,9 @@ class Instrument: ## Measured pseudo-ranging on TPS grid (high-frequency) logger.info("Computing measured pseudo-ranges on TPS") - self.tps_mprs = \ - self.local_timer_deviations \ - - self.distant_timer_deviations + self.ranging_noises + self.tps_mprs = self.the_wrt_tps_local - self.the_wrt_tps_distant \ + + self.ranging_noises + ## Test-mass interferometer local beams @@ -1203,12 +1205,12 @@ class Instrument: ## Sampling beatnotes, DWS measurements, and measured pseudo-ranges to THE grid - logger.info("Inverting timer deviations") - self.inverse_timer_deviations = self.local_timer_deviations \ - .transformed(lambda sc, x: self.invert_timer_deviations(x, sc)) + logger.info("Inverting THE with respect to TPS") + self.tps_wrt_the = self.the_wrt_tps_local \ + .transformed(lambda sc, x: self.invert_the_wrt_tps(x, sc)) self.timestamped = \ - lambda mosa, x: self.interpolate(x, -self.inverse_timer_deviations.for_each_mosa()[mosa]) + lambda mosa, x: self.interpolate(x, -self.tps_wrt_the.for_each_mosa()[mosa]) logger.info("Sampling inter-spacecraft beatnotes to THE grid") @@ -1413,7 +1415,7 @@ class Instrument: ## Closing simulation logger.info("Simulation complete") - def simulate_noises(self): + def simulate_noises(self): # pylint: disable=too-many-statements """Generate noise time series.""" ## Laser noise @@ -1453,6 +1455,23 @@ class Instrument: lambda _, x: x if np.isscalar(x) else x[self.initial_telemetry_physics_size:] ) + logger.debug("Integrating clock noise fluctuations") + + # Include initial telemetry time period + self.integrated_clock_noise_fluctuations_withinitial = \ + ForEachSC(lambda sc: + cumulative_trapezoid(np.broadcast_to( + self.clock_noise_fluctuations_withinitial[sc], + self.physics_size + self.initial_telemetry_physics_size), + dx=self.physics_dt, initial=0) + ) + + # Slice to only select physics period + self.integrated_clock_noise_fluctuations = \ + self.integrated_clock_noise_fluctuations_withinitial.transformed( + lambda _, x: x if np.isscalar(x) else x[self.initial_telemetry_physics_size:] + ) + ## Modulation noise logger.info("Generating modulation noise") @@ -1529,10 +1548,13 @@ class Instrument: noises.dws(self.physics_fs, self.physics_size, self.dws_asds[mosa]) ) - ## TCB synchronization noise + ## MOC time correlation noise - self.tcb_sync_noises = ForEachSC(lambda sc: - noises.tcb_sync(self.telemetry_fs, self.telemetry_size, self.sync_asds[sc]) + self.moc_time_correlation_noises = ForEachSC(lambda sc: + noises.moc_time_correlation( + self.telemetry_fs, + self.telemetry_size, + self.moc_time_correlation_asds[sc]) ) ## Angular jitters @@ -1719,37 +1741,37 @@ class Instrument: if not just_locked: raise RuntimeError(f"cannot apply locking conditions to remaining lasers '{list(dependencies.keys())}'") - def invert_timer_deviations(self, timer_deviations, sc): - """Invert timer deviations of a given spacecraft. + def invert_the_wrt_tps(self, the_wrt_tps, sc): + """Invert THE with respect to TPS of a given spacecraft. We recursively solve the implicit equation dtau(tau) = dtau_hat(tau - dtau(tau)) until the convergence criteria (tolerance) is met, or we exceed the maximum number of iterations. Args: - timer_deviations: array of timer deviations + the_wrt_tps: array of THEs with respect to TPS sc: spacecraft index """ - logger.debug("Inverting timer deviations for spacecraft %s", sc) + logger.debug("Inverting THE with respect to TPS for spacecraft %s", sc) logger.debug("Solving iteratively (tolerance=%s s, maxiter=%s)", self.clockinv_tolerance, self.clockinv_maxiter) # Drop samples at the edges to compute error - edge = min(100, len(timer_deviations) // 2 - 1) + edge = min(100, len(the_wrt_tps) // 2 - 1) error = 0 niter = 0 - next_inverse = timer_deviations + next_inverse = the_wrt_tps while not niter or error > self.clockinv_tolerance: if niter >= self.clockinv_maxiter: logger.warning("Maximum number of iterations '%s' reached for SC %s (error=%.2E)", niter, sc, error) break logger.debug("Starting iteration #%s", niter) inverse = next_inverse - next_inverse = self.interpolate(timer_deviations, -inverse) + next_inverse = self.interpolate(the_wrt_tps, -inverse) error = np.max(np.abs((inverse - next_inverse)[edge:-edge])) logger.debug("End of iteration %s, with an error of %.2E s", niter, error) niter += 1 - logger.debug("End of timer deviation inversion after %s iterations with an error of %.2E s", niter, error) + logger.debug("End of THE with respect to TCB inversion after %s iterations with an error of %.2E s", niter, error) return inverse def _write_attr(self, hdf5, *names): @@ -1847,6 +1869,9 @@ class Instrument: self.pprs.write(hdf5, 'pprs') self.d_pprs.write(hdf5, 'd_pprs') + logger.debug("Writing TPS with respect to TCB to '%s'", output) + self.tps_wrt_tcb.write(hdf5, 'tps_wrt_tcb') + logger.debug("Writing gravitational-wave responses to '%s'", output) self.gws.write(hdf5, 'gws') @@ -1863,7 +1888,12 @@ class Instrument: logger.debug("Writing clock noise to '%s'", output) self.clock_noise_offsets.write(hdf5, 'clock_noise_offsets') self.clock_noise_fluctuations.write(hdf5, 'clock_noise_fluctuations') - self.clock_noise_fluctuations_withinitial.write(hdf5, 'clock_noise_fluctuations_withinitial') + self.clock_noise_fluctuations_withinitial.write( + hdf5, 'clock_noise_fluctuations_withinitial') + self.integrated_clock_noise_fluctuations_withinitial.write( + hdf5, 'integrated_clock_noise_fluctuations_withinitial') + self.integrated_clock_noise_fluctuations.write( + hdf5, 'integrated_clock_noise_fluctuations') logger.debug("Writing modulation noise to '%s'", output) self.modulation_noises.write(hdf5, 'modulation_noises') @@ -1885,8 +1915,8 @@ class Instrument: self.oms_rfi_carrier_noises.write(hdf5, 'oms_rfi_carrier_noises') self.oms_rfi_usb_noises.write(hdf5, 'oms_rfi_usb_noises') - logger.debug("Writing TCB synchronization noise to '%s'", output) - self.tcb_sync_noises.write(hdf5, 'tcb_sync_noises') + logger.debug("Writing MOC time correlation noise to '%s'", output) + self.moc_time_correlation_noises.write(hdf5, 'moc_time_correlation_noises') logger.debug("Writing spacecraft angular jitter to '%s'", output) self.sc_jitter_phis.write(hdf5, 'sc_jitter_phis') @@ -1906,9 +1936,14 @@ class Instrument: self.local_usb_offsets.write(hdf5, 'local_usb_offsets') self.local_usb_fluctuations.write(hdf5, 'local_usb_fluctuations') - logger.debug("Writing local timer deviations to '%s'", output) - self.local_timer_deviations.write(hdf5, 'local_timer_deviations') - self.local_timer_deviations_withinitial.write(hdf5, 'local_timer_deviations_withinitial') + logger.debug("Writing local THE with respect to TPS to '%s'", output) + self.the_wrt_tps_local.write(hdf5, 'the_wrt_tps_local') + + logger.debug("Writing THE with respect to TCB to '%s'", output) + self.the_wrt_tcb_withinitial.write(hdf5, 'the_wrt_tcb_withinitial') + + logger.debug("Writing MOC time correlations to '%s'", output) + self.moc_time_correlations.write(hdf5, 'moc_time_correlations') logger.debug("Writing tilt-to-length couplings to '%s'", output) self.local_ttls.write(hdf5, 'local_ttls') @@ -1920,8 +1955,8 @@ class Instrument: self.distant_usb_offsets.write(hdf5, 'distant_usb_offsets') self.distant_usb_fluctuations.write(hdf5, 'distant_usb_fluctuations') - logger.debug("Writing propagated timer deviations to '%s'", output) - self.distant_timer_deviations.write(hdf5, 'distant_timer_deviations') + logger.debug("Writing propagated THEs with respect to TCB to '%s'", output) + self.the_wrt_tps_distant.write(hdf5, 'the_wrt_tcb_distant') logger.debug("Writing propagated adjacent beams to '%s'", output) self.adjacent_carrier_offsets.write(hdf5, 'adjacent_carrier_offsets') @@ -1990,8 +2025,8 @@ class Instrument: self.tps_rfi_usb_offsets.write(hdf5, 'tps_rfi_usb_offsets') self.tps_rfi_usb_fluctuations.write(hdf5, 'tps_rfi_usb_fluctuations') - logger.debug("Writing inverted timer deviations to '%s'", output) - self.inverse_timer_deviations.write(hdf5, 'inverse_timer_deviations') + logger.debug("Writing TPS with respect to THE to '%s'", output) + self.tps_wrt_the.write(hdf5, 'tps_wrt_the') logger.debug("Writing inter-spacecraft beatnotes sampled to THE grid to '%s'", output) self.the_isi_carrier_offsets.write(hdf5, 'the_isi_carrier_offsets') @@ -2095,9 +2130,6 @@ class Instrument: self.rfi_carriers.write(hdf5, 'rfi_carriers') self.rfi_usbs.write(hdf5, 'rfi_usbs') - logger.debug("Write timer deviations from TCB to '%s'", output) - self.tcb_timer_deviations.write(hdf5, 'tcb_timer_deviations') - logger.info("Closing measurement file '%s'", output) def plot_fluctuations(self, output=None, skip=0): diff --git a/lisainstrument/noises.py b/lisainstrument/noises.py index 57b854b0028152abff22b216c3f79f8ea2043f0b..a49ff2ee7d7b55e730c690930693a1997ab558c7 100644 --- a/lisainstrument/noises.py +++ b/lisainstrument/noises.py @@ -321,19 +321,20 @@ def dws(fs, size, asd): logger.debug("Generating DWS measurement (fs=%s Hz, size=%s, asd=%s rad/sqrt(Hz))", fs, size, asd) return violet(fs, size, 2 * pi * asd) -def tcb_sync(fs, size, asd): - """TCB synchronization noise. +def moc_time_correlation(fs, size, asd): + """MOC time correlation noise. - High-level noise model for the uncertainty we have in estimating - the TCB timer deviations, i.e., the equivalent TCB times for the + High-level noise model for the uncertainty we have in computing the MOC + time correlation (or time couples), i.e., the equivalent TCB times for the equally-sampled TPS timestamps. Assumed to be a white noise in timing, - S_tcbsnc(f) [s] = asd^2. + S_moc(f) [s] = asd^2. Args: asd: amplitude spectral density [s/sqrt(Hz)] """ - logger.debug("Generating TCB synchronization noise (fs=%s Hz, size=%s, asd=%s s/sqrt(Hz))", fs, size, asd) + logger.debug("Generating MOC time correlation noise (fs=%s Hz, size=%s, " + "asd=%s s/sqrt(Hz))", fs, size, asd) return white(fs, size, asd)