From e0f8ff0e91e71a82eade081724915383f1f3fa58 Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 09:59:49 -0600 Subject: [PATCH 01/10] Update documentation and reference for cluster_emcee_walkers function --- ssapy/utils.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index f7eaccd..8e25c7d 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -210,9 +210,11 @@ def cluster_emcee_walkers( chain, lnprob, lnprior, thresh_multiplier=1, verbose=False ): """ - Down-select emcee walkers to those with the largest mean posteriors + Down-select emcee walkers to those with the largest posterior mean. - Follows the algorithm of Hou, Goodman, Hogg et al. (2012) + Follows the algorithm of Hou, Goodman, Hogg et al. (2012), An affine-invariant + sampler for exoplanet fitting and discovery in radial velocity data. + The Astrophysical Journal, 745(2), 198. """ import emcee from distutils.version import LooseVersion From 2cf8d7cfa38f85f91ceaebdaebe2af736103d396 Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 10:00:40 -0600 Subject: [PATCH 02/10] Fix typo in return type --- ssapy/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index 8e25c7d..29b3d46 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -382,7 +382,7 @@ def regularize_default(particles, weights, num_particles_out=None, dimension=6): :type dimension: int, optional :return: Deltas from original particles and their weights - :rtype: (numpy.ndarray, numpy.ndarry) + :rtype: (numpy.ndarray, numpy.ndarray) """ num_particles_in, dim_in = particles.shape if num_particles_out is None: From dabfe3c02f7cd2343bc78b15b467ad6a261ee713 Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 10:24:39 -0600 Subject: [PATCH 03/10] Clarify docstring comment about units of function arguments --- ssapy/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index 29b3d46..6feaa63 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -770,7 +770,7 @@ def xyz_to_lb(x, y, z): def lb_to_tan(lb, b, mul=None, mub=None, lcen=None, bcen=None): """Convert lb-like coordinates & proper motions to orthographic tangent plane. - Everything is in radians. If mul is None (default), transformed + All units are in radians. If mul is None (default), transformed proper motions will not be returned. The tangent plane is always chosen so that +Y is towards b = 90, and +X is towards +lb. From 62e441c6874e5b1b7774926167e94973d6a823c0 Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 10:39:05 -0600 Subject: [PATCH 04/10] Clarify that covar is covariance matrix in docstring --- ssapy/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index 6feaa63..72819f9 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -963,8 +963,8 @@ def sigma_points(f, x, C, scale=1, fixed_dimensions=None): def unscented_transform_mean_covar(f, x, C, scale=1): - """Compute mean and covar using unscented transform given a transformation - f, a point x, and a covariance C. + """Compute mean and covariance matrix using unscented transform given a + transformation f, a point x, and a covariance C. This uses the sigma point convention from sigma_points. It assumes that f(sigma_points)[i] is f evaluated at the ith sigma point. If f does From ba3735bdaecdf1d4fe321525df30757c56a31835 Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 11:10:58 -0600 Subject: [PATCH 05/10] Add docstring to function converting GPS time to terrestrial time --- ssapy/utils.py | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index 72819f9..3ec20f7 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -990,10 +990,26 @@ def unscented_transform_mean_covar(f, x, C, scale=1): def _gpsToTT(t): - # Check if t is astropy Time object if so convert to GPS seconds if not assume GPS seconds. Convert to TT seconds by adding 51.184. - # Divide by 86400 to get TT days. - # Then add the TT time of the GPS epoch, expressed as an MJD, which - # is 44244.0 + """ + Convert GPS time to Terrestrial Time (TT). + + Parameters: + ---------- + t : Time or float + If `t` is an instance of `astropy.time.Time`, it is assumed to represent GPS time in the form of an Astropy Time object. + If `t` is a float, it is assumed to be GPS time in seconds. + + Returns: + ------- + float + The equivalent time in Terrestrial Time (TT) expressed in days since the GPS epoch (January 6, 1980). + + Notes: + ----- + The conversion involves adding a constant offset of 51.184 seconds to the GPS time, + then converting the result into days by dividing by 86400 (the number of seconds in a day). + The epoch for GPS is represented as a Modified Julian Date (MJD) of 44244.0. + """ if isinstance(t, Time): t = t.gps return 44244.0 + (t + 51.184) / 86400 From 9189426cbd13b79feb0c9e4b38c779d57dc12ddd Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 12:56:57 -0600 Subject: [PATCH 06/10] Remove redundant comment --- ssapy/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index 3ec20f7..41ab470 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -1187,7 +1187,7 @@ def unit_vector(vector): return vector / np.linalg.norm(vector) -def get_angle(a, b, c): # a,b,c where b is the vertex +def get_angle(a, b, c): """ Calculate the angle between two vectors where b is the vertex of the angle. From 2b774e8d8eaf6a6214cfd55a0b284efbbca34f4f Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 13:04:24 -0600 Subject: [PATCH 07/10] Remove second definition of function that appears twice in utils.py --- ssapy/utils.py | 39 --------------------------------------- 1 file changed, 39 deletions(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index 41ab470..c6bdbed 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -2129,42 +2129,3 @@ def find_smallest_bounding_cube(r): upper_bound = center + half_side_length return lower_bound, upper_bound - - -def points_on_circle(r, v, rad, num_points=4): - # Convert inputs to numpy arrays - r = np.array(r) - v = np.array(v) - - # Find the perpendicular vectors to the given vector v - if np.all(v[:2] == 0): - if np.all(v[2] == 0): - raise ValueError("The given vector v must not be the zero vector.") - else: - u = np.array([1, 0, 0]) - else: - u = np.array([-v[1], v[0], 0]) - u = u / np.linalg.norm(u) - w = np.cross(u, v) - w_norm = np.linalg.norm(w) - if w_norm < 1e-15: - # v is parallel to z-axis - w = np.array([0, 1, 0]) - else: - w = w / w_norm - # Generate a sequence of angles for equally spaced points - angles = np.linspace(0, 2 * np.pi, num_points, endpoint=False) - - # Compute the x, y, z coordinates of each point on the circle - x = rad * np.cos(angles) * u[0] + rad * np.sin(angles) * w[0] - y = rad * np.cos(angles) * u[1] + rad * np.sin(angles) * w[1] - z = rad * np.cos(angles) * u[2] + rad * np.sin(angles) * w[2] - - # Apply rotation about z-axis by 90 degrees - rot_matrix = np.array([[0, 1, 0], [-1, 0, 0], [0, 0, 1]]) - rotated_points = np.dot(rot_matrix, np.column_stack((x, y, z)).T).T - - # Translate the rotated points to the center point r - points_rotated = rotated_points + r.reshape(1, 3) - - return points_rotated \ No newline at end of file From af781d33ebdf06d52655feb2579d14bfaca8024c Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 13:22:33 -0600 Subject: [PATCH 08/10] Rename conversion functions with _to_ convention --- ssapy/utils.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index c6bdbed..f62024b 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -1538,7 +1538,7 @@ def lonlat_distance(lat1, lat2, lon1, lon2): return (c * EARTH_RADIUS) -def altitude2zenithangle(altitude, deg=True): +def altitude_to_zenithangle(altitude, deg=True): if deg: out = 90 - altitude else: @@ -1546,7 +1546,7 @@ def altitude2zenithangle(altitude, deg=True): return out -def zenithangle2altitude(zenith_angle, deg=True): +def zenithangle_to_altitude(zenith_angle, deg=True): if deg: out = 90 - zenith_angle else: @@ -1554,7 +1554,7 @@ def zenithangle2altitude(zenith_angle, deg=True): return out -def rightasension2hourangle(right_ascension, local_time): +def rightascension_to_hourangle(right_ascension, local_time): if type(right_ascension) is not str: right_ascension = dd_to_hms(right_ascension) if type(local_time) is not str: @@ -1570,7 +1570,7 @@ def rightasension2hourangle(right_ascension, local_time): def equatorial_to_horizontal(observer_latitude, declination, right_ascension=None, hour_angle=None, local_time=None, hms=False): if right_ascension is not None: - hour_angle = rightasension2hourangle(right_ascension, local_time) + hour_angle = rightascension_to_hourangle(right_ascension, local_time) if hms: hour_angle = hms_to_dd(hour_angle) elif hour_angle is not None: @@ -1587,7 +1587,7 @@ def equatorial_to_horizontal(observer_latitude, declination, right_ascension=Non zenith_angle = np.arccos(np.sin(observer_latitude) * np.sin(declination) + np.cos(observer_latitude) * np.cos(declination) * np.cos(hour_angle)) - altitude = zenithangle2altitude(zenith_angle, deg=False) + altitude = zenithangle_to_altitude(zenith_angle, deg=False) _num = np.sin(declination) - np.sin(observer_latitude) * np.cos(zenith_angle) _den = np.cos(observer_latitude) * np.sin(zenith_angle) @@ -1602,7 +1602,7 @@ def equatorial_to_horizontal(observer_latitude, declination, right_ascension=Non def horizontal_to_equatorial(observer_latitude, azimuth, altitude): altitude, azimuth, latitude = np.radians([altitude, azimuth, observer_latitude]) - zenith_angle = zenithangle2altitude(altitude) + zenith_angle = zenithangle_to_altitude(altitude) zenith_angle = [-zenith_angle if latitude < 0 else zenith_angle][0] From a6b25556dd7d417dc430e2bf4a0985196972ff7c Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 13:38:08 -0600 Subject: [PATCH 09/10] Add docstring to lonlat_distance function --- ssapy/utils.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/ssapy/utils.py b/ssapy/utils.py index f62024b..2efb078 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -1528,14 +1528,32 @@ def ra_dec(r=None, v=None, x=None, y=None, z=None, vx=None, vy=None, vz=None, r_ def lonlat_distance(lat1, lat2, lon1, lon2): - # Haversine formula + """Calculate the great-circle distance between two points + on Earth's surface using the Haversine formula. + + Parameters + ---------- + lat1 : float + Latitude of the first point in radians. + lat2 : float + Latitude of the second point in radians. + lon1 : float + Longitude of the first point in radians. + lon2 : float + Longitude of the second point in radians. + + Returns + ------- + distance : float + Distance between the two points in kilometers. + """ dlon = lon2 - lon1 dlat = lat2 - lat1 a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2 c = 2 * np.arcsin(np.sqrt(a)) - # Radius of earth in kilometers. Use 3956 for miles - # calculate the result - return (c * EARTH_RADIUS) + # Calculate distance in kilometers, use 3956 for miles + distance = c * EARTH_RADIUS + return distance def altitude_to_zenithangle(altitude, deg=True): From 0b76bc6704fefce6710097ab9d456c0605125800 Mon Sep 17 00:00:00 2001 From: Jason Bernstein Date: Wed, 25 Sep 2024 13:54:28 -0600 Subject: [PATCH 10/10] Add newline between Parameters and Returns sections of docstring --- ssapy/utils.py | 1 + 1 file changed, 1 insertion(+) diff --git a/ssapy/utils.py b/ssapy/utils.py index 2efb078..7884244 100644 --- a/ssapy/utils.py +++ b/ssapy/utils.py @@ -2028,6 +2028,7 @@ def get_times(duration, freq, t): The starting time. Default is "2025-01-01". example input: duration=(30, 'day'), freq=(1, 'hr'), t=Time("2025-01-01", scale='utc') + Returns ------- times: array-like