From 8fcae94da77b9618f54f19ad118df4c53d438dc1 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 28 Aug 2026 10:53:48 +0200 Subject: [PATCH 1/2] Add SpEC's OmegaDotEccRemoval.py --- docs/Api.rst | 5 + pyproject.toml | 1 + .../EccentricityControlParams.py | 181 ++ .../EccentricityControl/OmegaDotEccRemoval.py | 1960 +++++++++++++++++ .../Test_EccentricityControlParams.py | 59 + 5 files changed, 2206 insertions(+) create mode 100644 src/SimulationSupport/EccentricityControl/EccentricityControlParams.py create mode 100644 src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py create mode 100644 tests/EccentricityControl/Test_EccentricityControlParams.py diff --git a/docs/Api.rst b/docs/Api.rst index c78e672..31a2807 100644 --- a/docs/Api.rst +++ b/docs/Api.rst @@ -5,3 +5,8 @@ API Reference :members: :undoc-members: :show-inheritance: + +.. automodule:: SimulationSupport.EccentricityControl.EccentricityControlParams + :members: + :undoc-members: + :show-inheritance: diff --git a/pyproject.toml b/pyproject.toml index de6be2c..7cfd303 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,6 +18,7 @@ dependencies = [ "gpytorch", "pandas", "sxs", + "varpro @ git+https://github.com/sxs-collaboration/varpro.git@978106eaf3d7a6a7f0c5f167726d8e0fc59fc95d", ] [project.optional-dependencies] diff --git a/src/SimulationSupport/EccentricityControl/EccentricityControlParams.py b/src/SimulationSupport/EccentricityControl/EccentricityControlParams.py new file mode 100644 index 0000000..c68561c --- /dev/null +++ b/src/SimulationSupport/EccentricityControl/EccentricityControlParams.py @@ -0,0 +1,181 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. +"""Estimate eccentricity and updated orbital parameters from trajectories.""" + +import logging +from pathlib import Path +from typing import Dict, Literal, Optional, Union + +import numpy as np + +from .OmegaDotEccRemoval import ( + ComputeOmegaAndDerivsFromFile, + FindTmin, + performAllFits, +) + +logger = logging.getLogger(__name__) + +# Keys of the dictionary returned by 'eccentricity_control_params' +EccentricityParams = Literal[ + "Eccentricity", + "EccentricityError", + "Omega0", + "Adot0", + "D0", + "DeltaOmega0", + "DeltaAdot0", + "DeltaD0", + "NewOmega0", + "NewAdot0", + "NewD0", + "Tmin", + "Tmax", +] + + +def eccentricity_control_params( + trajectory_a: np.ndarray, + trajectory_b: np.ndarray, + separation: float, + orbital_angular_velocity: float, + radial_expansion_velocity: float, + mass_a: float, + mass_b: float, + spin_a: Optional[np.ndarray] = None, + spin_b: Optional[np.ndarray] = None, + tmin: Optional[float] = None, + tmax: Optional[float] = None, + target_eccentricity: float = 0.0, + plot_output_dir: Optional[Union[str, Path]] = None, +) -> Dict[EccentricityParams, float]: + r"""Get new orbital parameters for a binary system to control eccentricity. + + The eccentricity is estimated from the trajectories of the binary objects + and updates to the orbital parameters are suggested to drive the orbit to + the target eccentricity, using SpEC's ``OmegaDotEccRemoval.py``. Currently + supports only circular target orbits (target eccentricity = 0). + + Parameters + ---------- + trajectory_a : numpy.ndarray + Trajectory of the first object, with shape ``(num_times, 4)``. The + first column is the time, the remaining three columns are the + coordinates. + trajectory_b : numpy.ndarray + Trajectory of the second object, in the same format as + ``trajectory_a``. + separation : float + Initial coordinate separation ``D_0`` of the two objects, i.e. the + initial data parameter that is being controlled. + orbital_angular_velocity : float + Initial orbital angular velocity ``Omega_0``. + radial_expansion_velocity : float + Initial radial expansion velocity ``adot_0``. + mass_a : float + Christodoulou mass of the first object at the reference time. + mass_b : float + Christodoulou mass of the second object at the reference time. + spin_a : numpy.ndarray, optional + Dimensionful spin of the first object, with shape ``(num_times, 4)``. + The first column is the time, the remaining three columns are the + components of the spin vector. If either spin is unspecified, spin + effects are ignored in the fits. + spin_b : numpy.ndarray, optional + Dimensionful spin of the second object, in the same format as + ``spin_a``. + tmin : float, optional + The lower time bound for the eccentricity estimate. Used to remove + initial junk and transients in the data. If unspecified, uses SpEC's + ``OmegaDotEccRemoval.FindTmin`` to estimate it. + tmax : float, optional + The upper time bound for the eccentricity estimate. A reasonable value + would include 2-3 orbits. Default is ``500 + 5 * pi / Omega_0``. + target_eccentricity : float, optional + Eccentricity that the updated orbital parameters should achieve. + Currently only 0 (circular orbit) is supported. + plot_output_dir : str or pathlib.Path, optional + Output directory for plots. + + Returns + ------- + dict + Dictionary with the keys listed in ``EccentricityParams``. + """ + if target_eccentricity != 0.0: + raise ValueError( + "Only circular orbits are currently supported for eccentricity" + " control." + ) + + # Compute the orbital frequency and its time derivative from the + # trajectories + t, Omega, dOmegadt, OmegaVec = ComputeOmegaAndDerivsFromFile( + trajectory_a, trajectory_b + ) + + # Set time bounds if not provided + if tmin is None: + tmin = max(FindTmin(t, dOmegadt, 500), t[0]) + if tmax is None: + tmax = min(500 + 5 * np.pi / orbital_angular_velocity, t[-1]) + logger.info( + "Estimating eccentricity from trajectory data in time range" + f" {tmin:.3f} to {tmax:.3f}." + ) + + # Call into SpEC's OmegaDotEccRemoval.py + eccentricity, delta_Omega0, delta_adot0, delta_D0, ecc_std_dev, _ = ( + performAllFits( + XA=trajectory_a, + XB=trajectory_b, + t=t, + Omega=Omega, + dOmegadt=dOmegadt, + OmegaVec=OmegaVec, + mA=mass_a, + mB=mass_b, + sA=spin_a, + sB=spin_b, + IDparam_omega0=orbital_angular_velocity, + IDparam_adot0=radial_expansion_velocity, + IDparam_D0=separation, + tmin=tmin, + tmax=tmax, + tref=tmin, + opt_freq_filter=True, + opt_varpro=True, + opt_type="bbh", + opt_tmin=tmin, + opt_improved_Omega0_update=True, + check_periastron_advance=True, + plot_output_dir=plot_output_dir, + Source="", + ) + ) + logger.info( + f"Eccentricity estimate is {eccentricity:g} +/- {ecc_std_dev:e}." + " Update orbital parameters as follows" + f" for target eccentricity {target_eccentricity:g} (choose two):\n" + f"Omega0 += {delta_Omega0:e} ->" + f" {orbital_angular_velocity + delta_Omega0:.8g}\n" + f"adot0 += {delta_adot0:e} ->" + f" {radial_expansion_velocity + delta_adot0:e}\n" + f"D0 += {delta_D0:e} -> {separation + delta_D0:.8g}" + ) + # These keys must correspond to 'EccentricityParams' + return { + "Eccentricity": eccentricity, + "EccentricityError": ecc_std_dev, + "Omega0": orbital_angular_velocity, + "Adot0": radial_expansion_velocity, + "D0": separation, + "DeltaOmega0": delta_Omega0, + "DeltaAdot0": delta_adot0, + "DeltaD0": delta_D0, + "NewOmega0": orbital_angular_velocity + delta_Omega0, + "NewAdot0": radial_expansion_velocity + delta_adot0, + "NewD0": separation + delta_D0, + "Tmin": tmin, + "Tmax": tmax, + } diff --git a/src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py b/src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py new file mode 100644 index 0000000..5f819fa --- /dev/null +++ b/src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py @@ -0,0 +1,1960 @@ +# fmt: off +# isort: skip_file + +""" +This code was copied from the SpEC implementation +('Support/DatDataManip/OmegaDotEccRemoval.py') with only minimal changes and +should be modernized. It was kept in its original form to make transitioning +from SpEC to this package easier. + +Performs eccentricity fitting, and outputs the estimated eccentricity +and updated initial data parameters. +""" + +from __future__ import division +import matplotlib +import sys, os +import matplotlib.pyplot as plt +import numpy as np +from numpy import sqrt, size, sin, cos, array, pi, cross, dot, mean +from numpy.random import rand +import argparse +import re +from scipy import optimize, signal +from scipy.interpolate import CubicSpline, InterpolatedUnivariateSpline +import varpro as VarPro +from varpro import RankError +from io import StringIO + + +########################################################################## +### Functions copied from other SpEC modules +### +### 'Support/Python/Utils.py': ReadH5, ReadDat, error, warning, norm, +### and SmoothData +### 'Support/Python/ParseIdParams.py': ParseIdParams and ParseIdParamsDict +### 'Support/Python/BbhDiagnosticsImpl.py': Compute_OrbitalFrequency +########################################################################## + + +def ReadH5(file): + """Opens the file 'file' as read-only H5 file, and returns the +h5py.File object. +""" + import h5py + try: + f = h5py.File(file, 'r') + except IOError: + os.sys.stderr.write("\n############# ERROR #############\n") + os.sys.stderr.write('Attempt to open h5-file "' + file + '" failed.\n') + raise + return f + + +################################################################ + + +def ReadDat(file, cols=None): + """Read an ASCII .dat file, and return the specified columns. +file -- filename of the file to read +cols -- columns to read, in the format expected from numpy.loadtxt's usecols + option, namely an iterable listing the columns, like (1,4,6,5) + If ==None, return all columns +RETURNS: + a numpy array +""" + from numpy import loadtxt + try: + return loadtxt(file, usecols=cols, ndmin=2) # always return 2D + except: + os.sys.stderr.write("\n############# ERROR #############\n") + os.sys.stderr.write('Attempt to read "' + file + + '" as dat file failed.\n') + raise + + +################################################################ + + +def error(msg): + e_msg = "\n############# ERROR #############\n{}\n".format(msg) + raise Exception(e_msg) + + +################################################################ + + +# Print a warning message to stderr +def warning(msg): + os.sys.stderr.write("\n############# WARNING #############\n") + os.sys.stderr.write("{}\n".format(msg)) + + +################################################################ + + +def norm(data, axis=None): + """Returns the L2-norm over a specific axis. The axis option +behaves exactly like in numpy.sum: if axis=None (default), then +sum over all axes. If negative, count from the last axis. +""" + import numpy as np + return np.sqrt(np.sum(np.square(data), axis=axis)) + + +################################################################ + + +def SmoothData(t, + data, + width, + Deriv=0, + conv_factor=4., + cut_off=None, + Tstart=None, + DT=None, + dt_tol=1e-12): + """Smooth the data (t,data) using Gaussian convolution. +Derivatives do not seem to work well when timesteps are much larger than 1. +It also does not work (though it should) for unequally spaced timesteps. +The smooth data will be truncated on both sides based on the smoothing width. +""" + import math + import numpy as np + + # Convert to expected objects + width = float(width) + try: + data.shape + except AttributeError: + data = np.array(data) + if len(data.shape) == 1: + data = data[:, np.newaxis] #promote the array with singleton dim + if len(data.shape) != 2: + error("data is expected to be 2d array, not {}d".format(len( + data.shape))) + + try: + t.shape + except AttributeError: + t = np.array(t) + + convolution_width = width * conv_factor + if cut_off == None: + cut_width = width * conv_factor + else: + cut_width = width * cut_off + + # compute average time-spacing + def ComputeAvgDt(t, dt_tol): + dt = np.diff(t) + aveDt = np.mean(dt) + MaxDtDiff = max(abs(dt - aveDt)) + if MaxDtDiff > dt_tol: + error( + "Times are not equally spaced within dt_tol={}.\n".format( + dt_tol) + + "Maximum deviation from the average is {}.".format(MaxDtDiff)) + return aveDt + + aveDt = ComputeAvgDt(t, dt_tol) + + # set up moving window + if (t[-1] - t[0]) <= 6 * width: + error("Total time-interval must be at least 6x larger than width") + + # setup output times + if DT == None: + t_out = t + else: + Tstart = t[0] if (Tstart is None) else Tstart + if Tstart < t[0]: + warn("--Tstart smaller than first time-step %g" % t[0]) + t_out = np.arange(Tstart, t[-1], step=DT) + + def K(x, width): + return 1. / np.sqrt(math.pi) / width * np.exp(-x * x / (width * width)) + + def K1(x, width): + return -2. * x / np.sqrt(np.pi) / width**3 * np.exp(-x * x / + (width * width)) + + def K2(x, width): + return (4*x*x-2*width**2)/np.sqrt(np.pi)/width**5 \ + * np.exp(-x*x/(width**2)) + + idxmin = 0 + idxmax = 0 + + # limit t_out to range so that the kernel fits + t_out = t_out[(t_out >= t[0] + cut_width) * (t_out <= t[-1] - cut_width)] + + output_array = np.empty((len(t_out), 1 + len(data[0]))) + for curr_idx in range(len(t_out)): + curr_t = t_out[curr_idx] + # update bounds for convolution + while (t[idxmin] < curr_t - convolution_width): + idxmin = idxmin + 1 + while idxmax < len(t) - 1 and t[idxmax] < curr_t + convolution_width: + idxmax = idxmax + 1 + + if (t[idxmax] > curr_t + convolution_width): + idxmax = idxmax - 1 + + if (abs(curr_t - t[0]) < cut_width or abs(curr_t - t[-1]) < cut_width): + continue + + if (idxmax == idxmin): + error("For t=%f, no data-points in interval.\n" % curr_t) + + # compute sums + # s = Sum Kf, w=Sum K + f12 = K(curr_t - t[idxmin:idxmax + 1], width) + w = np.sum(aveDt * (f12[0:-1] + f12[1:])) + s = np.sum(aveDt * (data[idxmin:idxmax] * f12[0:-1].reshape( + (-1, 1)) + data[idxmin + 1:idxmax + 1] * f12[1:].reshape((-1, 1))), + axis=0) + + # s1= Sum K1f, w1=Sum K1 + if (Deriv >= 1): + f12 = K1(curr_t - t[idxmin:idxmax + 1], width) + w1 = np.sum(aveDt * (f12[0:-1] + f12[1:])) + s1 = np.sum(aveDt * (data[idxmin:idxmax] * f12[0:-1].reshape( + (-1, 1)) + data[idxmin + 1:idxmax + 1] * f12[1:].reshape( + (-1, 1))), + axis=0) + + # s2= Sum K2f, w2=Sum K2 + if (Deriv >= 2): + f12 = K2(curr_t - t[idxmin:idxmax + 1], width) + w2 = np.sum(aveDt * (f12[0:-1] + f12[1:])) + s2 = np.sum(aveDt * (data[idxmin:idxmax] * f12[0:-1].reshape( + (-1, 1)) + data[idxmin + 1:idxmax + 1] * f12[1:].reshape( + (-1, 1))), + axis=0) + + ### combine for answer ### + if (Deriv == 0): + result = np.concatenate(((curr_t, ), s / w)) + elif (Deriv == 1): + result = np.concatenate(((curr_t, ), s1 / w - s / w**2 * w1)) + elif (Deriv == 2): + # Note: This formula produces visible noise + # for non-uniformly spaced data even when the + # 'w1' terms are included, which sould take + # correct for non-uniform spacing. + # (Harald, Nov 29, 2006) + result = np.concatenate(((curr_t, ), s2 / w - 2. * s1 * w1 / w**2 - + s * w2 / w**2 + 2. * s * w1 * w1 / w**3)) + output_array[curr_idx, :] = result[:] + return output_array + + +################################################################ + + +def ParseIdParams(path,file="ID_Params.perl",array=False): + """Parses an ID_Params.perl file and returns a dictionary of {key: value} +pairs. Arguments are: + path # path to file to parse + file=ID_Params.perl # filename of file to parse + array=False # if True, output values in array context, not string + # context, where + # string context => 'value1, value 2' + # array context => ['value1','value2'] +""" + + fullpath = os.path.join(path,file) + try: + text = open(fullpath, 'r').read() + except IOError: + error("ParseIdParams: Could not read file %s!" %fullpath) + + BaseRegexp = r'([^\s]+)\s*= *([^\n]*)' + Regexp = r"^(?:@|\$)%s;\s*$" %BaseRegexp + + # Grab all (key,value) pairs and put them in a dictionary + PairList = re.findall(Regexp, text, re.MULTILINE) + PairDict = dict(PairList) + for key in PairDict: + # strip leading/trailing spaces, overall parens, and ALL quotes + PairDict[key] = PairDict[key].strip(' ()').replace('"','').replace("'",'') + if array: + for key in PairDict: + PairDict[key] = [i.strip() for i in PairDict[key].split(',')] + + return ParseIdParamsDict(PairDict) + +# Returns a more useful error message when dict key does not exist +class ParseIdParamsDict(dict): + def __getitem__(self, key): + try: + val = dict.__getitem__(self, key) + except KeyError: + error("Key '{}' does not exist in this ParseIdParams dict.\n" + "Available keys are {}".format(key,self.keys())) + return val + + +################################################################ + + +def Compute_OrbitalFrequency(xA, xB, N, method="Fit", NSamples=None): + r"""Given numpy arrays for the black hole locations xA, xB, + perform fits to 2N+1 data-points around each point, and from + the fit compute the instantaneous orbital frequency, + Omega = r\times \dot{r} / r^2 + return t,Omega + """ + + def FitData(data, N, NSamples=None): + """given a numpy array data with time as first column, perform fits + covering N points before and after each data-point for each column. + return the fitted values and their first time-derivatives as a + numpy arrays, with first column time""" + # collect output data, first column being time + last_idx=len(data)-N-1 + + if NSamples==None: + step=1 + else: + step=max(int(last_idx/NSamples),1) + + # The output times + t_final = data[N:last_idx:step,0] + + x_tmp = [] + v_tmp = [] + for idx in range(N, last_idx, step): + # Time at which we want the result + T = data[idx,0] + + x = data[idx-N:idx+N+1,0]-T # Shift back to t=0 + y = data[idx-N:idx+N+1,1:] # Fit all the columns at once! + p0 = np.polyfit(x,y,2) + x_tmp.append(p0[2]) # p0[2] is the constant part of fit + v_tmp.append(p0[1]) # p0[1] is the linear part + + return np.column_stack((t_final,x_tmp)), np.column_stack((t_final,v_tmp)) + + def FilterData(data,N): + assert False, "The FilterData function is poorly tested" + from scipy.ndimage import gaussian_filter1d + mode = 'reflect' + return gaussian_filter1d(data,N,mode=mode),\ + gaussian_filter1d(data,N,mode=mode,order=1) + + def SmoothData(data,N): + from SimulationSupport.EccentricityControl.OmegaDotEccRemoval import ( + SmoothData, + ) + t = data[:,0] + dat = data[:,1:] + dt_interp = None + return SmoothData(t,dat,N,DT=dt_interp), \ + SmoothData(t,dat,N,DT=dt_interp,Deriv=1) + + def SplineData(data): + """Interpolating spline (no smoothing)""" + from scipy.interpolate import splrep,splev + t = data[:,0] + spline_x = splrep(t, data[:,1]) + spline_y = splrep(t, data[:,2]) + spline_z = splrep(t, data[:,3]) + dx = splev(t, spline_x, der=1) + dy = splev(t, spline_y, der=1) + dz = splev(t, spline_z, der=1) + v = np.vstack((t,dx,dy,dz)).T + return data, v + + if NSamples is not None and method!='Fit': + error("NSamples only works with 'Fit'") + + if method=="Fit": + data = np.column_stack((xA,xB[:,1:4])) + xs_fit,vs = FitData(data,N,NSamples=NSamples) + xA_fit = xs_fit[:,0:4] + xB_fit = xs_fit[:,[0,4,5,6]] + vA = vs[:,0:4] + vB = vs[:,[0,4,5,6]] + elif method=="Filter": + xA_fit,vA=FilterData(xA,N) + xB_fit,vB=FilterData(xB,N) + elif method=="Smooth": + xA_fit,vA=SmoothData(xA,N) + xB_fit,vB=SmoothData(xB,N) + elif method=="Spline": + xA_fit,vA=SplineData(xA) + xB_fit,vB=SplineData(xB) + else: + error("Don't know method '{}'".format(method)) + + # Compute Orbital frequency (r x dr/dt)/r^2 + t = xA_fit[:,0] + dr = xA_fit[:,1:] - xB_fit[:,1:] + dr2 = norm(dr,axis=1)**2 #slightly inefficient + dv = vA[:,1:] - vB[:,1:] + Omega = np.cross(dr,dv)/dr2[:,np.newaxis] + return t,Omega + + +########################################################################## +### Functions for unloading input and calculating related quantities +########################################################################## + +# Choose the earliest tmin **after** the junk radiation. +# This function computes the absolute value of the second derivative of +# omega, ODblDot; it then takes a running average of this abs. value +# over avgp points; separately, it averages ODblDot from time 500 until the +# end and calls this yfin; with this, it calculates drop which is the value of +# ODblDot @t=0 - yfin; a threshold of ODblDot, ythr, is the sum of yfin and +# threshp*drop; then tmin is the time at which ODblDot first goes under this +# threshold, ythr, + tshiftp(set @200); however, tmin is capped at 500; +# in other words, tmin=min(tmin, 500). (Author: Robert McGehee) +def FindTmin(t_temp2, dOmegadt, max_tmin): + # Compute the 2nd deriv. of Omega, take its abs. value + d2Omegadt = (dOmegadt[2:]-dOmegadt[0:-2])/(t_temp2[2:]-t_temp2[0:-2]) + t=t_temp2[1:-1] + ODblDot=abs(d2Omegadt) + + # Compute the running average of ODblDot with a given avg param + avgp=10 #How many points in the running average (has to be an even #) + y= [0.00]*size(t) + for i in range(0,avgp//2): + for j in range(0,avgp+1): + if j!=i: + y[i]=y[i]+ODblDot[j] + y[i]=y[i]/avgp + for i in range(avgp//2,size(t)-avgp//2): + for j in range(i-avgp//2,i+avgp//2+1): + if j!=i: + y[i]=y[i]+ODblDot[j] + y[i]=y[i]/avgp + for i in range(size(t)-avgp//2,size(t)): + for j in range(size(t)-avgp-1,size(t)): + if j!=i: + y[i]=y[i]+ODblDot[j] + y[i]=y[i]/avgp + + # Compute the average,yfin, of ODblDot from t=500 until the end + sum=0.0000000 + for i in range(500,size(t)): + sum=sum+ODblDot[i] + yfin=sum/(size(t)-500) + drop=ODblDot[0]-yfin #the drop in magnitude determines a threshold + threshp=0.001 #parameter which controls how low the threshold ythr is + ythr=yfin+threshp*drop + + # Loop through y to find first time where y < ythr + k=0 + while (k=ythr): + k=k+1 + if k==size(t): + k=k-1 + tshiftp=200 #parameter controls how much "safety" is added to tmin + tmin=t[k]+tshiftp + + # Don't allow tmin to be arbitrarily large + return min(tmin, max_tmin) + +#==== logic for idperl option +# Get Omega0, adot0, D0, as well as the initial ADM masses +# of neutron stars. Return nan for BH masses (we should calculate +# them using GetRelaxedMasses instead) +def ParseIDParams(path, binary_type): + File = os.path.basename(os.path.realpath(path)) + Dir = os.path.dirname(os.path.realpath(path)) + D = ParseIdParams(Dir, file=File, array=True) + + Omega0 = float(D['ID_Omega0'][0]) + adot0 = float(D['ID_adot0'][0]) + D0 = float(D['ID_d'][0]) + + mA = float("nan") + mB = float("nan") + + if binary_type=="bhns": + mB = float(D['ADMmassNS'][0]) + if binary_type=="nsns": + mA = float(D['ADMmassNS1'][0]) + mB = float(D['ADMmassNS2'][0]) + + return Omega0,adot0,D0,mA,mB + +def GetTrajectories(Dir, Type, tmax_fit): + if Type == "bbh": + Horizons = ReadH5(os.path.join(Dir,"Horizons.h5")) + TrajA = Horizons["AhA.dir/CoordCenterInertial.dat"] + TrajB = Horizons["AhB.dir/CoordCenterInertial.dat"] + elif Type == "bhns": + Horizons = ReadH5(os.path.join(Dir,"Horizons.h5")) + Matter = ReadH5(os.path.join(Dir,"Matter.h5")) + TrajA = Horizons["AhA.dir/CoordCenterInertial.dat"] + TrajB = Matter["InertialCenterOfMassNS1.dat"] + elif Type == "nsns": + Matter = ReadH5(os.path.join(Dir,"Matter.h5")) + TrajA = Matter["InertialCenterOfMassNS1.dat"] + TrajB = Matter["InertialCenterOfMassNS2.dat"] + + # Truncate the trajectories to a region around the fit interval + tmax = 1.5 * tmax_fit + tA = TrajA[:,0] + tB = TrajB[:,0] + # Some versions of h5py (version 3.0, 3.1 ?) cannot be indexed with + # boolean arrays. https://github.com/h5py/h5py/issues/1847 + # So here we create integer arrays instead of saying + # things like TrajA[tA= 0.6*IDparam_omega0 and + f <= 1.4*IDparam_omega0 for f in peak_freqs]] + + # Pick peak of FFT in frequency space to be guess for omega_0. Only consider + # frequencies within 40% of the initial data omega0. If there is more than + # one peak or there are no peaks in this range, default to 0.8*omega0 as a + # guess and don't filter. + if len(viable_peaks) != 1: + return 0.8*IDparam_omega0, None, None + # Filtering can only be done if a minimum occurs after the guess peak + # frequency. To filter we set the frequency spectrum for frequencies higher + # than the first minimum after guess peak to 0. + elif len(minimum_freqs[minimum_freqs > viable_peaks[0]]) > 0: + filter_threshold = minimum_freqs[minimum_freqs > viable_peaks[0]][0] + filtered_FFT = np.fft.rfft(pre_FFT) + filtered_FFT[FFT_freq > filter_threshold] = 0 + filtered_dOmegadt = (np.fft.irfft(filtered_FFT, n=signal_length) + .real[pad_length:-pad_length] / window) \ + + detrend_diff + return viable_peaks[0], t_grid, filtered_dOmegadt + else: + return viable_peaks[0], None, None + +########################################################################## +### Functions for performing fits and computing updates +########################################################################## + +# fit the data (t,y) to the model F[p,t], by least-squares fitting the params +# p. +def fit(t, y, F, p0, bounds, jac, name): + + # t,y -- arrays of equal length, having t and y values of data + # F -- function: F(p,t) taking parameters p and set of t-values + # p0 -- starting values of parameters + residual = lambda p,t,y: F(p,t) - y + # the jacobian function passed to least_squares must have the form + # jac(p, t, y) + jacfunc = lambda p,t,y: [jac(p,t_i) for t_i in t] + res = optimize.least_squares( + residual, + method = "dogbox", + x0 = p0, + jac = jacfunc, + bounds = bounds, + args = (t,y), + ) + + if not res.success: + if name == "F2cos2_SS": + error("minimize failed in {}: {}".format(name, res.message)) + else: + warning("minimize failed in {}: {}".format(name, res.message)) + + return res.x, res.cost, res.success + +#TODO: this should be updated now that we bound omega +def CheckPeriastronAdvance(omega,Omega0,name,ecc,ecc_str,summary): + # See discussion just before section IV of arXiv:1012.1549. The + # idea here is that omega/Omega0 should be exactly equal to 1 for a + # pure Newtonian orbit, but for GR or PN orbits omega/Omega0 should + # be slightly less than 1 because of periastron advance. + # Furthermore, omega/Omega0 should be independent of eccentricity + # for small eccentricity. + # + # So if omega/Omega0 is much larger than 1, then this is unphysical + # because it represents negative periastron advance. + # If omega/Omega0 is much smaller than 1, then this is also very large + # periastron advance and unphysical (one can argue where the cutoff should + # go). + # + # The big question is what should we do if we encounter these + # unphysical situations. Currently we flag for a human to look at + # by setting ecc=9.99999. Note that such a large value of eccentricity + # will currently cause EccReduce.pm to report a failure. + # + # These situations might occur when + # eccentricity is so small that we cannot find an accurate fit. + # Currently we have no reliable way of detecting inaccurate fits + # (other than the res/B check); but with varpro we will have an + # error estimate for ecc, so perhaps the periastron check will be + # encountered less often in that case. + # + if omega/Omega0<0.5 or omega/Omega0>1.2: + ecc = 9.99999 + ecc_str = str(ecc) + err="OmegaDotEccRemoval: {name}-fit resulted in large (>1.2) or " \ + "negative (<0.5) periastron advance:\nomega/Omega0={omegafrac}." \ + " This is likely wrong (or a bad omega fit from too small ecc " \ + "oscillation amplitudes),\nso eccentricity was set to {ecc}\n" \ + .format(name=name,omegafrac=omega/Omega0,ecc=ecc) + warning(err) + summary.write(err) + return ecc,ecc_str + +def ComputeUpdate(Omega0, adot0, D0, + Tc, B, omega, + phi, phi_tref, + name, tmin, tmax, rms, + Improved_Omega0_update, check_periastron_advance, + params_output_dir, Source, summary, + B_std_dev = None, omega_std_dev = None): + ''' + Computes and returns eccentricity and corrections to be added to initial + data parameters Omega_0, adot, and D. Corrected initial data removes + spurious eccentricity from the evolution. Only two out of the three + corrections returned should be applied to initial data before restarting + evolution. See arXiv:1012.1549 for derivation. + ''' + + delta_adot0=B/(2.*Omega0)*cos(phi) + delta_Omega0=-B*omega/(4.*Omega0*Omega0)*sin(phi) + if(Improved_Omega0_update): + # extra factor Omega0/omega in delta_Omega0 + delta_Omega0=-B/(4.*Omega0)*sin(phi) + + delta_D0=-B*D0*omega*sin(phi)/(2*Omega0*(Omega0**2+2/D0**3)) + ecc=B/(2.*Omega0*omega) + ecc_str="{:7.7f}".format(ecc) + ecc_std_dev = None + if (B_std_dev != None) and (omega_std_dev != None): + ecc_std_dev = sqrt(B**2 * omega_std_dev**2 + omega**2 * B_std_dev**2) \ + / (2 * Omega0 * omega * omega) + + if rms/B>0.4: + # See discussion just before section IV of arXiv:1012.1549. + # The idea here is that if the oscillations are very small + # (meaning that B is small enough that it is smaller than the residual + # of the fit), then the fit effectively does not see the oscillations. + # So the eccentricity estimate is not very accurate, and we report + # a bound here. + # + # When we switch to varpro, we will presumably have an error bound + # on B, and instead of rms/B > 0.4 we can use some criterion like + # (error bound on B) < B. + summary.write("Large residual of ecc-fit, report bound on ecc\n") + # for a sine-wave, the rms is 1/2 its amplitude. Therefore, assume + # that the amplitude of a bad-fit is 2*rms. Double that, for safety + # and because we have to disregard the term omega/Omega0 + ecc=4.*rms/(2.*Omega0*Omega0) + ecc_str="<{:.1e}".format(ecc) + elif check_periastron_advance: + ecc,ecc_str = CheckPeriastronAdvance(omega,Omega0,name, + ecc,ecc_str,summary) + pi=np.pi + summary.write("%s: %+11.8f %+11.8f " \ + "%+9.6f %s %9.6f %9.6f\n" + %(name,delta_Omega0,delta_adot0,delta_D0,ecc_str, + (phi-pi/2.)%(2.*pi), + (phi_tref-pi/2.)%(2.*pi) # mean anomaly, in [0, 2pi] + )) + + if params_output_dir: + f = open(os.path.join(params_output_dir, "Params_%s.dat"%name), 'w') + f.write("# EccRemoval.py utilizing orbital frequency "\ + "Omega, fit %s\n"%name) + f.write("# Source file=%s\n" % Source) + f.write("# Omega0=%10.8g, adot0=%10.8g, D0=%10.8g\n"%(Omega0,adot0,D0)) + f.write("# Fitting interval [tmin,tmax]=[%g,%g]\n"%(tmin,tmax)) + f.write("# oscillatory part of fit: (B,omega,phi)=(%g,%g,%g)\n" + %(B,omega,phi)) + f.write("# ImprovedOmega0Update=%s\n"%Improved_Omega0_update) + f.write("# [1] = Omega0\n") + f.write("# [2] = 1e4 adot0\n") + f.write("# [3] = D0\n") + f.write("# [4] = ecc\n") + f.write("%10.12f\t%10.12f\t%10.12f\t%10.12f\n" + %(Omega0, 1e4*adot0, D0, ecc)) + f.write("%10.12f\t%10.12f\t%10.12f\t%10.12f\n" + %(Omega0+delta_Omega0, 1e4*(adot0+delta_adot0), + D0+delta_D0, 0.)) + f.close() + + f = open("Fit_%s.dat"%name, 'w') + f.write("# EccRemoval.py utilizing orbital frequency Omega, fit %s\n" + %name) + f.write("# Source file=%s\n" % Source) + f.write("# Omega0=%10.8g, adot0=%10.8g\n"%(Omega0,adot0)) + f.write("# Fitting interval [tmin,tmax]=[%g,%g]\n"%(tmin,tmax)) + f.write("# oscillatory part of fit: (B,omega,phi)=(%g,%g,%g)\n" + %(B,omega,phi)) + f.write("# [1] = Tstart\n") + f.write("# [2] = Tend\n") + f.write("# [3] = Tc\n") + f.write("# [4] = B\n") + f.write("# [5] = omega\n") + f.write("# [6] = sin(phi)\n") + f.write("# [7] = rms residual\n") + f.write("# [8] = rms residual/B\n") + f.write("# [9] = omega/Omega0\n") + f.write("# [10] = ecc\n") + + f.write("%g %g %g %g %g %g %g %g %g %g\n" + %(tmin, tmax, Tc, B, omega, sin(phi), + rms, rms/B, omega/Omega0, ecc)) + f.close() + + return ecc, delta_Omega0, delta_adot0, delta_D0, ecc_std_dev + +def plot_fitted_function(p, F, x, y, idxFit, idxPlot, idxZoom, name, style, + axes): + # add the plot to all four panels + # p, F: params and fitting function + # x, y - complete data + # idxFit -- indices used in fit, for indicating fit-interval + # idxPlot -- indices to be plotted in left windows + # idxZoom -- indices to be plotted in right windows + # name -- string placed into legend + # style - plot-style + + xBdry=array([x[idxFit][0], x[idxFit][-1]]) + yBdry=array([y[idxFit][0], y[idxFit][-1]]) + + ((ax1,ax2),(ax3,ax4)) = axes + + ## Top-left plot: fit + ax1.plot(x[idxPlot],1e6*F(p,x[idxPlot]),style, label=name) + + # add point at begin and end of fit interval: + ax1.plot(xBdry,1e6*F(p,xBdry),'o') + + # bottom-left plot: residual + data=1e6*(F(p,x)-y) + ax3.plot(x[idxPlot],data[idxPlot],style,label=name) + ax3.plot(xBdry,1e6*(F(p,xBdry)-yBdry),'o') + ylim3=ax3.get_ylim() + miny=min(ylim3[0], min(data[idxFit])) + maxy=max(ylim3[1], max(data[idxFit])) + ax3.set_ylim(miny, maxy) + ax3.set_title("1e6 residual") + + ## Top-right plot: zoom of fit + ax2.plot(x[idxZoom],1e6*F(p,x[idxZoom]),style, label=name) + ax2.legend(loc='upper right', bbox_to_anchor=(1.15,1.3) + ,labelspacing=0.25,handletextpad=0.0,fancybox=True + ) + + # add point at begin of fit interval: + ax2.plot([xBdry[0]],[1e6*F(p,xBdry)[0]],'o') + + # bottom-right plot: zoom of residual + ax4.plot(x[idxZoom],1e6*(F(p,x[idxZoom])-y[idxZoom]),style,label=name) + ax4.plot([xBdry[0]],[(1e6*(F(p,xBdry)-yBdry))[0]],'o') + ax4.set_title("1e6 residual") + + return + +def make_full_plot(Source,t,dOmegadt,tmin,tmax,idxFit, + params,funcs,names,linestyles,plotfilename): + #=== set plotting intervals + idxPlot= (t> tmin - 0.2*(tmax-tmin)) & (t< tmax + 0.35*(tmax-tmin) ) + idxZoom=(t< tmin + 0.2*(tmax-tmin)) + tPlot=t[idxPlot] + dOmegadtPlot=dOmegadt[idxPlot] + + fig,axes = plt.subplots(2,2) + fig.text(0.5, 0.95, + "%s [%g, %g]" % (os.path.basename(Source), tmin, tmax), + color='b', size='large', ha='center') + + ((ax1,ax2),(ax3,ax4)) = axes + + #=== Top left plot dOmega/dt === + ax1.plot(tPlot[tPlot >= tmin], 1e6 * dOmegadtPlot[tPlot >= tmin], 'k', + label="dOmega/dt", linewidth=2) + xlim1=ax1.get_xlim() + ylim1=ax1.get_ylim() + ax1.plot(tPlot, 1e6 * dOmegadtPlot, 'k',label="dOmega/dt", linewidth=2) + xlim2=ax1.get_xlim() + ax1.set_xlim(xlim2) + ax1.set_ylim(ylim1) + ax1.set_title("1e6 dOmega/dt") + + #==== bottom left ==== + # set x-axes to top-left scale, and y-axes to something small + # as initial conditions for adding line by line above + ax3.set_xlim(xlim1) + ax3.set_ylim([-1e-10, 1e-10]) + + #==== Top right plot -- zoom of dOmega/dt ==== + ax2.plot(t[idxZoom], 1e6 * dOmegadt[idxZoom], + 'k', linewidth=2) # label="dOmega/dt", + ax2.set_title("1e6 dOmega/dt .") # extra space to avoid legend + + # Plot individual fits + for (p,f,name,linestyle) in zip(params,funcs,names,linestyles): + plot_fitted_function(p, f, t, dOmegadt, idxFit, idxPlot, + idxZoom, name, linestyle, axes) + + # zoom out of the y-axis in the lower left panel by 15% + ylim1 = ax3.get_ylim() + Deltay=ylim1[1]-ylim1[0] + ax3.set_ylim([ylim1[0]-0.15*Deltay, ylim1[1]+0.15*Deltay]) + + # adjust margins of figure + fig.subplots_adjust(left=0.09, right=0.95, bottom=0.07, hspace=0.25) + + plt.savefig(plotfilename) + +class FitBounds: + """Set bounds for some of the variables being fit""" + def __init__(self, tmax, IDparam_omega0): + # Time to coalescence must be greater than the fit interval + self.Tc = [tmax+2.0*np.spacing(abs(tmax)), np.inf] + # Frequency should be positive and similar to the initial frequency + self.omega = [0.6*IDparam_omega0, 1.4*IDparam_omega0] + # Amplitude of the cos term is chosen to be positive (negative is phi+=pi) + self.B = [1e-16, np.inf] + # phi is a phase, trivially bounded + self.phi = [0, 2*pi] + # For various variables that have no limits + self.none = [-np.inf, np.inf] + +def performNonspinFits(t,dOmegadt,idxFit,tFit,dOmegadtFit, + tmin,tmax,tref,IDparam_omega0, + IDparam_D0,IDparam_adot0,q,omega_guess,opt_tmin, + opt_improved_Omega0_update,check_periastron_advance, + params_output_dir,plot_output_dir,Source,summary): + ''' + Fit the eccentricity estimator given in arXiv:1012.1549 to the given + trajectory, not including spin terms, and compute initial data from the + resulting best fit parameters. A series of fits is performed, each one + including a new term from the estimator model. Returns eccentricity and + initial data corrections based on the best fit to the complete model. + ''' + + # Set bounds for some of the variables + lim = FitBounds(tmax, omega_guess) + + # 0PN approximation + Tmerger = 5. / (64.*(1. - ((q-1.) / (q+1.))**2.)*IDparam_omega0**(8./3.)) + + # fit a0(T-t)^(-11./8) + F1 = lambda p,t: p[1]*(p[0]-t)**(-11/8) + p0 = [Tmerger, 1e-5] + pBounds = list(zip(*[lim.Tc, lim.none])) + jac = lambda p,t: [ (-11/8)*p[1]*(p[0]-t)**(-19/8), + (p[0]-t)**(-11/8) ] + pF1, rmsF1, F1_status = fit(tFit, dOmegadtFit, F1, p0, pBounds, jac, "F1") + + #================ + + F1cos1 = lambda p,t: p[1]*(p[0]-t)**(-11/8) + p[2]*cos(p[3]*t+p[4]) + jac = lambda p,t: [ (-11/8)*p[1]*(p[0]-t)**(-19/8), + (p[0]-t)**(-11/8), + cos(p[3]*t + p[4]), + -p[2]*sin(p[3]*t + p[4])*t, + -p[2]*sin(p[3]*t + p[4]) ] + + rmsF1cos1 = 2*rmsF1 + pF1cos1 = pF1 + # try a few initial guesses for the phase to ensure good convergence + for phi in range(0,6): + p0 = [pF1[0], pF1[1], rmsF1, omega_guess, phi] + pBounds = list(zip(*[lim.Tc, lim.none, lim.B, lim.omega, lim.phi])) + ptemp, rmstemp, F1cos1_status = fit(tFit, dOmegadtFit, F1cos1, + p0, pBounds, jac, "F1cos1") + if(rmstemp t[-1]: + # After this point, tmax can no longer be arbitrarily large. + # tmax must be inside the dataset. + tmax = t[-1] + + # Get masses of compact objects + # (At tmin for BHs, at t=0 for NSs) + mA, mB = GetRelaxedMasses(Source, tmin, opts.t, mA_id, mB_id) + + if opts.t=="bbh": + Horizons = ReadH5(os.path.join(opts.d, "Horizons.h5")) + sA = Horizons["AhA.dir/DimensionfulInertialSpin.dat"] + sB = Horizons["AhB.dir/DimensionfulInertialSpin.dat"] + else: + sA = None + sB = None + + params_output_dir = None if opts.no_output_files else "." + plot_output_dir = None if opts.no_plot else "." + + _, _, _, _, _, summary_string = performAllFits(IDparam_omega0, IDparam_adot0, + IDparam_D0, XA, XB, mA, mB, sA, + sB, t, Omega, dOmegadt, OmegaVec, + tmin, tmax, tref, + opts.freq_filter, opts.varpro, + opts.t, opts.tmin, + opts.improved_Omega0_update, + not opts.no_check, + params_output_dir, + plot_output_dir, Source, + True) + + if not opts.no_output_files: + summary_file = open("summary.txt", "w") + summary_file.write(summary_string) + summary_file.close() + + print(summary_string) + +def performAllFits(IDparam_omega0, IDparam_adot0, IDparam_D0, XA, XB, mA, mB, + sA, sB, t, Omega, dOmegadt, OmegaVec, tmin, tmax, tref, + opt_freq_filter, opt_varpro, opt_type, opt_tmin, + opt_improved_Omega0_update, check_periastron_advance, + params_output_dir=None, plot_output_dir=None, Source=None, + always_do_spin_fits=False): + ''' + Fit an eccentricity estimator to the given trajectory. Return the + eccentricity, final initial data corrections, and output for printing. + + If variable projection is enabled, also return estimated uncertainty of + eccentricity. Otherwise, return None. + ''' + + summary = StringIO() + q = max(mA/mB, mB/mA) # q > 1 by convention + + idxFit= (t>=tmin) & (t<=tmax) + tFit=t[idxFit] + OmegaFit=Omega[idxFit] + dOmegadtFit=dOmegadt[idxFit] + + omega_guess, filtered_t, filtered_dOmegadt = \ + computeOmegaGuessAndFilterTraj(tFit, dOmegadtFit, IDparam_omega0) + if (opt_freq_filter and filtered_t is not None): + tFit = filtered_t + dOmegadtFit = filtered_dOmegadt + + # First, do the usual fits + if opt_varpro: + nonspin_ecc, nonspin_delta_Omega0, nonspin_delta_adot0, \ + nonspin_delta_D0, nonspin_ecc_std_dev = \ + performNonspinVarPro(t,dOmegadt,idxFit,tFit,dOmegadtFit,tmin,tmax,tref, + IDparam_omega0,IDparam_D0,IDparam_adot0,q, + omega_guess,opt_tmin,opt_improved_Omega0_update, + check_periastron_advance,params_output_dir, + plot_output_dir,Source,summary) + else: + nonspin_ecc, nonspin_delta_Omega0, nonspin_delta_adot0, \ + nonspin_delta_D0, nonspin_ecc_std_dev = \ + performNonspinFits(t,dOmegadt,idxFit,tFit,dOmegadtFit,tmin,tmax,tref, + IDparam_omega0,IDparam_D0,IDparam_adot0,q,omega_guess, + opt_tmin,opt_improved_Omega0_update, + check_periastron_advance,params_output_dir, + plot_output_dir,Source,summary) + + if ((opt_freq_filter and filtered_t is not None + and not always_do_spin_fits) + or not (opt_type=="bbh" and sA is not None and sB is not None)): + # If filtering is turned on and has succeeded, spin contribution is removed + # from trajectory before fitting, so return corrections from nonspin fits + # is not necessary. So return early. + # Also return early if we aren't going to do spin fits anyway because + # we are not BBH and we don't have spins. + summary_string = summary.getvalue() + summary.close() + return nonspin_ecc, nonspin_delta_Omega0, nonspin_delta_adot0, \ + nonspin_delta_D0, nonspin_ecc_std_dev, summary_string + + # Now, do fits that include spin-spin interactions, + # see Buonnano et al., 2010 (arXiv 1012.1549v2). + + # Everything below is only for BBH (since NS have negligible spins) + alpha_intrp, S_0_perp_n, T_merge, Amp = \ + GetVarsFromSpinData(sA, sB, XA, XB, mA, mB, OmegaVec, t, tmin) + + if opt_varpro: + spin_ecc, spin_delta_Omega0, spin_delta_adot0, spin_delta_D0, \ + spin_ecc_std_dev = \ + performSpinVarPro(t,dOmegadt,idxFit,tFit,dOmegadtFit,alpha_intrp,tmin, + tmax,tref,IDparam_omega0,IDparam_D0,IDparam_adot0, + q,omega_guess,opt_tmin,opt_improved_Omega0_update, + check_periastron_advance,params_output_dir, + plot_output_dir,Source,summary) + else: + spin_ecc, spin_delta_Omega0, spin_delta_adot0, spin_delta_D0, \ + spin_ecc_std_dev = \ + performSpinFits(t,dOmegadt,idxFit,tFit,OmegaFit,dOmegadtFit,T_merge, + Amp,alpha_intrp,S_0_perp_n,tmin,tmax,tref, + IDparam_omega0,IDparam_D0,IDparam_adot0,omega_guess, + opt_tmin,opt_improved_Omega0_update, + check_periastron_advance,params_output_dir, + plot_output_dir,Source,summary) + + summary_string = summary.getvalue() + summary.close() + + return spin_ecc, spin_delta_Omega0, spin_delta_adot0, spin_delta_D0, \ + spin_ecc_std_dev, summary_string + +#============================================================================= + +if __name__ == "__main__": + + # Use a non-interactive backend when run as a script + matplotlib.use('Agg') + + p = argparse.ArgumentParser( + description=__doc__, + formatter_class=argparse.RawTextHelpFormatter + ) + + required = p.add_argument_group("required arguments") + required.add_argument("--idperl", type=str, required=True, metavar="FILE", + help="File with initial data variables ID_{Omega0,adot0,D0},\n"\ + "e.g. ID_Params.perl") + + required.add_argument("-d", type=str, required=True, metavar="DIR", + help="Directory containing the evolution trajectory file(s)") + required.add_argument("-t", type=str, required=True, + choices=['bbh','bhns','nsns'], + help="""Type of binary determines which files are expected in DIR. +bbh: + Horizons.h5/Ah{A,B}.dir/CoordCenterInertial.dat + Horizons.h5/Ah{A,B}.dir/DimensionfulInertialSpin.dat + Horizons.h5/Ah{A,B}.dir/ChristodoulouMass.dat +bhns: + Horizons.h5/AhA.dir/CoordCenterInertial.dat + Matter.h5/InertialCenterOfMassNS1.dat +nsns: + Matter.h5/InertialCenterOfMassNS{1,2}.dat +""") + + p.add_argument("--tmin",type=float, metavar="FLOAT", + help="""Fit points with t>tmin. +Default determined by FindTmin estimation scheme.""") + p.add_argument("--tmax",type=float, metavar="FLOAT", + help="""Fit points with t Date: Wed, 2 Sep 2026 09:53:00 +0200 Subject: [PATCH 2/2] Fix a missing square in OmegaDotEccRemoval.py --- src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py b/src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py index 5f819fa..59212a6 100644 --- a/src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py +++ b/src/SimulationSupport/EccentricityControl/OmegaDotEccRemoval.py @@ -1512,7 +1512,7 @@ def nonspinFPhiF0(alpha, t): # Calculate quantities for correction formulae B_nonspin = np.sqrt(pLin_nonspin[2]**2 + pLin_nonspin[3]**2) B_std_dev = sqrt(pLin_nonspin[2]**2 * param_std_dev_nonspin[2]**2 \ - + pLin_nonspin[3]**2 * param_std_dev_nonspin[3]) / B_nonspin + + pLin_nonspin[3]**2 * param_std_dev_nonspin[3]**2) / B_nonspin phi_nonspin = np.arctan2(pLin_nonspin[3], pLin_nonspin[2]) # phi should be constrained to (0, 2pi), # but np.arctan2 is constrained to (-pi, pi)