Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,4 @@ src/pygenray/__about__.py
docs/_build/
dist/
.vscode/
uv.lock
uv.lock
18 changes: 18 additions & 0 deletions .pre-commit-config.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
repos:
- repo: https://github.com/pre-commit/pre-commit-hooks
rev: "v6.0.0"
hooks:
- id: check-added-large-files
args: ["--maxkb=5120"]
- id: check-merge-conflict
- id: check-yaml
- id: end-of-file-fixer
- id: trailing-whitespace
- id: debug-statements

- repo: https://github.com/astral-sh/ruff-pre-commit
rev: "v0.15.11"
hooks:
- id: ruff-check
args: ["--fix", "--show-fixes"]
- id: ruff-format
1 change: 0 additions & 1 deletion .readthedocs.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -23,4 +23,3 @@ python:
path: .
extra_requirements:
- docs

2 changes: 1 addition & 1 deletion docs/building_notes.md
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,4 @@ actually... scipy might be worth it. Because I can have the ray code done in jus
- use launched ray and previous bracketing rays and repeat until convergence


**where I left off**: I have built out `shoot_rays` and `shoot_ray`, and created a single method that can be called from within the functions that take an environment as input, and when accessing shared memory. But I haven't tested all of that functionality yet.
**where I left off**: I have built out `shoot_rays` and `shoot_ray`, and created a single method that can be called from within the functions that take an environment as input, and when accessing shared memory. But I haven't tested all of that functionality yet.
2 changes: 1 addition & 1 deletion docs/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
import sys

# Add the source directory to the path so Sphinx can import the modules
sys.path.insert(0, os.path.abspath('../src'))
sys.path.insert(0, os.path.abspath("../src"))

project = "pygenray"
copyright = "2025, John Ragland"
Expand Down
1 change: 0 additions & 1 deletion docs/quick_start.md
Original file line number Diff line number Diff line change
@@ -1 +0,0 @@

2 changes: 1 addition & 1 deletion docs/ray_physics.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,4 +4,4 @@
- Environment depth and source / receiver depths are specified with positive z, defined by distance from surface of ocean in meters.
- The ray state, $y$, defined by $y = \left[\mathrm{T, z, p_z} \right]^T$, uses sign convention where z is negative and increases in z correspond to moving towards the surface.
- Subsequently, a positive ray angle corresponds to a ray moving towards the surface.
- A positive ray ID corresponds to a positive launch angle (towards the surface)
- A positive ray ID corresponds to a positive launch angle (towards the surface)
3 changes: 3 additions & 0 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,9 @@ exclude_lines = [
"if TYPE_CHECKING:",
]

[tool.ruff.lint.per-file-ignores]
"src/pygenray/__init__.py" = ["F403"]

[dependency-groups]
dev = [
"pytest>=8.3.5",
Expand Down
166 changes: 120 additions & 46 deletions src/pygenray/eigenrays.py
Original file line number Diff line number Diff line change
@@ -1,25 +1,27 @@
"""
Tools and methods for calculating eigenrays for specifed receiver depths.
"""

import pygenray as pr
import numpy as np
import multiprocessing as mp
from tqdm import tqdm


def find_eigenrays(
rays,
receiver_depths,
source_depth,
source_range,
receiver_range,
num_range_save,
environment,
ztol=1,
max_iter=20,
num_workers=None,
**kwargs,
):
'''
rays,
receiver_depths,
source_depth,
source_range,
receiver_range,
num_range_save,
environment,
ztol=1,
max_iter=20,
num_workers=None,
**kwargs,
):
"""
Given an initial ray fan, find eigenrays with [regula falsi](https://en.wikipedia.org/wiki/Regula_falsi#The_regula_falsi_(false_position)_method) method of root finding.

Parameters
Expand Down Expand Up @@ -51,7 +53,7 @@ def find_eigenrays(
-------
erays : dict
dictionary of eigen rays. Key values are values in `receiver_depths`.
'''
"""
erays_dict = {}
num_eigenrays = {}
num_eigenrays_found = {}
Expand All @@ -60,7 +62,7 @@ def find_eigenrays(
for rd_idx, receiver_depth in enumerate(receiver_depths):
## get initial bracketing rays
# get indices before sign changes
depth_sign = np.sign(rays.zs[:,-1] + receiver_depth)
depth_sign = np.sign(rays.zs[:, -1] + receiver_depth)
sign_change = np.diff(depth_sign)
bracket_idxs_start = np.where(sign_change)[0]

Expand All @@ -70,11 +72,11 @@ def find_eigenrays(
bracket_idxs = np.column_stack([bracket_idxs_start, bracket_idxs_start + 1])

# compute regula falsi launch angles
z1s = rays.zs[bracket_idxs[:,0].astype(int),-1]
z2s = rays.zs[bracket_idxs[:,1].astype(int),-1]
z1s = rays.zs[bracket_idxs[:, 0].astype(int), -1]
z2s = rays.zs[bracket_idxs[:, 1].astype(int), -1]

theta1s = rays.thetas[bracket_idxs[:,0].astype(int)]
theta2s = rays.thetas[bracket_idxs[:,1].astype(int)]
theta1s = rays.thetas[bracket_idxs[:, 0].astype(int)]
theta2s = rays.thetas[bracket_idxs[:, 1].astype(int)]

erays_dict[rd_idx] = []
failed_eray_theta_brackets[rd_idx] = []
Expand All @@ -90,7 +92,7 @@ def find_eigenrays(
for eray_idx in range(num_eigenrays[receiver_depth]):
z1_distance = np.abs(z1s[eray_idx] + receiver_depth)
z2_distance = np.abs(z2s[eray_idx] + receiver_depth)

if (z1_distance < ztol) or (z2_distance < ztol):
bracket_winner = np.argmin([z1_distance, z2_distance])

Expand All @@ -105,31 +107,55 @@ def find_eigenrays(
eray_found_idxs.append(eray_idx)
else:
continue

# remove found eigenrays from search
z1s = np.delete(z1s, eray_found_idxs)
z2s = np.delete(z2s, eray_found_idxs)
theta1s = np.delete(theta1s, eray_found_idxs)
theta2s = np.delete(theta2s, eray_found_idxs)
"""

regula_falsi_thetas = theta1s - (z1s + receiver_depth) * (theta2s - theta1s) / (z2s - z1s)
regula_falsi_thetas = theta1s - (z1s + receiver_depth) * (theta2s - theta1s) / (
z2s - z1s
)

if len(regula_falsi_thetas) > 5: # use parallel processing for large number of rays
if (
len(regula_falsi_thetas) > 5
): # use parallel processing for large number of rays
# construct argument iterable for parallel processing
args_list = []
for k in range (len(regula_falsi_thetas)):
args = (k, z1s[k], z2s[k], theta1s[k], theta2s[k], regula_falsi_thetas[k],
receiver_depth, source_depth, source_range, receiver_range,
num_range_save, environment, ztol, max_iter, kwargs)
for k in range(len(regula_falsi_thetas)):
args = (
k,
z1s[k],
z2s[k],
theta1s[k],
theta2s[k],
regula_falsi_thetas[k],
receiver_depth,
source_depth,
source_range,
receiver_range,
num_range_save,
environment,
ztol,
max_iter,
kwargs,
)
args_list.append(args)

# map individual eigen ray finding to different workers
if num_workers is None:
num_workers = mp.cpu_count()
with mp.get_context('spawn').Pool(num_workers) as pool:
results = list(tqdm(pool.imap(_find_single_eigenray, args_list), total=len(args_list), desc="Finding eigenrays"))

with mp.get_context("spawn").Pool(num_workers) as pool:
results = list(
tqdm(
pool.imap(_find_single_eigenray, args_list),
total=len(args_list),
desc="Finding eigenrays",
)
)

# Filter out None results and add successful rays
for result in results:
if result is not None:
Expand All @@ -138,10 +164,26 @@ def find_eigenrays(
failed_eray_theta_brackets[rd_idx].append((theta1s[k], theta2s[k]))

else: # use sequential processing for small number of rays
for k in tqdm(range(len(regula_falsi_thetas)), desc='Finding eigenrays:'):
ray = _find_single_eigenray((k, z1s[k], z2s[k], theta1s[k], theta2s[k], regula_falsi_thetas[k],
receiver_depth, source_depth, source_range, receiver_range,
num_range_save, environment, ztol, max_iter, kwargs))
for k in tqdm(range(len(regula_falsi_thetas)), desc="Finding eigenrays:"):
ray = _find_single_eigenray(
(
k,
z1s[k],
z2s[k],
theta1s[k],
theta2s[k],
regula_falsi_thetas[k],
receiver_depth,
source_depth,
source_range,
receiver_range,
num_range_save,
environment,
ztol,
max_iter,
kwargs,
)
)
if ray is not None:
erays_dict[rd_idx].append(ray)
else:
Expand All @@ -150,28 +192,59 @@ def find_eigenrays(
num_eigenrays_found[rd_idx] = len(erays_dict[rd_idx])

# Create EigenRays object after processing all receiver depths
erays = pr.EigenRays(receiver_depths, erays_dict, environment, num_eigenrays, num_eigenrays_found, failed_eray_theta_brackets)
erays = pr.EigenRays(
receiver_depths,
erays_dict,
environment,
num_eigenrays,
num_eigenrays_found,
failed_eray_theta_brackets,
)
return erays


def _find_single_eigenray(args):
"""
Find single Eigen ray given the bracketing ray depths, and launch angles.
"""
k, z1, z2, theta1, theta2, regula_falsi_theta, receiver_depth, source_depth, source_range, receiver_range, num_range_save, environment, ztol, max_iter, kwargs = args

(
k,
z1,
z2,
theta1,
theta2,
regula_falsi_theta,
receiver_depth,
source_depth,
source_range,
receiver_range,
num_range_save,
environment,
ztol,
max_iter,
kwargs,
) = args

iter_count = 0
# Regula Falsi root finding loop
while True:

ray = pr.shoot_ray(source_depth, source_range, regula_falsi_theta, receiver_range, num_range_save, environment, **kwargs)
ray = pr.shoot_ray(
source_depth,
source_range,
regula_falsi_theta,
receiver_range,
num_range_save,
environment,
**kwargs,
)

if ray is None:
print(f'Failed to find eigen ray for receiver depth {receiver_depth} [m] and approximate launch angle {regula_falsi_theta} [m] ray θ = 90°')
print(
f"Failed to find eigen ray for receiver depth {receiver_depth} [m] and approximate launch angle {regula_falsi_theta} [m] ray θ = 90°"
)
return None

if np.abs(ray.z[-1] + receiver_depth) < ztol:

if np.abs(ray.z[-1] + receiver_depth) < ztol:
# flip launch angle to match sign convention
ray.launch_angle = -ray.launch_angle
return ray
Expand All @@ -185,13 +258,14 @@ def _find_single_eigenray(args):
z2 = ray.z[-1]
theta2 = regula_falsi_theta

regula_falsi_theta = theta1 - (z1 + receiver_depth) * (theta2 - theta1) / (z2 - z1)
regula_falsi_theta = theta1 - (z1 + receiver_depth) * (theta2 - theta1) / (
z2 - z1
)

if iter_count > max_iter:
return None

iter_count += 1

iter_count += 1


__all__ = ['find_eigenrays']
__all__ = ["find_eigenrays"]
Loading
Loading