From d0d2c819900eb760cd9d7d11b98df33085409c58 Mon Sep 17 00:00:00 2001 From: John Ragland Date: Wed, 18 Mar 2026 15:42:51 -0400 Subject: [PATCH 1/2] added tests --- .github/workflows/tests.yml | 27 +++ .gitignore | 3 +- pyproject.toml | 6 + tests/conftest.py | 55 +++++ tests/fixtures/munk_regression.npz | Bin 0 -> 7586 bytes tests/test_environment.py | 180 ++++++++++++++++ tests/test_physics.py | 324 +++++++++++++++++++++++++++++ tests/test_ray_objects.py | 278 +++++++++++++++++++++++++ 8 files changed, 872 insertions(+), 1 deletion(-) create mode 100644 .github/workflows/tests.yml create mode 100644 tests/conftest.py create mode 100644 tests/fixtures/munk_regression.npz create mode 100644 tests/test_environment.py create mode 100644 tests/test_physics.py create mode 100644 tests/test_ray_objects.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..f49e72e --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,27 @@ +name: Tests + +on: + push: + branches: [main] + pull_request: + +jobs: + test: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + + - name: Install uv + uses: astral-sh/setup-uv@v5 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.12" + + - name: Install dependencies + run: uv sync --group dev + + - name: Run tests + run: uv run pytest tests/ diff --git a/.gitignore b/.gitignore index 1ca5add..e33c601 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ __pycache__/ .DS_Store docs/_build/ dist/ -.vscode/ \ No newline at end of file +.vscode/ +uv.lock \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 5a6cded..72c9be9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -74,3 +74,9 @@ exclude_lines = [ "if __name__ == .__main__.:", "if TYPE_CHECKING:", ] + +[dependency-groups] +dev = [ + "pytest>=8.3.5", + "pytest-cov>=5.0.0", +] diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..4c25277 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,55 @@ +""" +Shared test fixtures for pygenray tests. +""" +import matplotlib +matplotlib.use('Agg') + + +def pytest_addoption(parser): + """Register --regenerate-physics CLI flag for physics regression tests.""" + parser.addoption( + '--regenerate-physics', action='store_true', default=False, + help='Regenerate physics regression fixture and skip comparison.', + ) + +import numpy as np +import pytest +import xarray as xr + +from pygenray.ray_objects import Ray, RayFan + + +def _make_ray(launch_angle: float, source_depth: float, n_bottom: int = 0, + n_surface: int = 0, N: int = 10, R: float = 10000.0) -> Ray: + """Helper: build a synthetic Ray without running the ODE solver. + + The y array uses the positive-z convention expected by Ray.__init__: + y[0,:] = travel times + y[1,:] = depth (positive = deeper) + y[2,:] = ray parameter sin(θ)/c (positive for downward ray in ODE) + """ + r = np.linspace(0.0, R, N) + t = r / 1500.0 + # Depth increases linearly (simulating a shallow downward ray) + z_ode = np.linspace(source_depth, source_depth + R * 0.01, N) + p_ode = np.ones(N) * np.sin(np.radians(abs(launch_angle) + 1e-3)) / 1500.0 + y = np.vstack([t, z_ode, p_ode]) # shape (3, N) + return Ray(r=r, y=y, n_bottom=n_bottom, n_surface=n_surface, + launch_angle=launch_angle, source_depth=source_depth) + + +@pytest.fixture +def simple_ray(): + """Single synthetic Ray with 10 range steps.""" + return _make_ray(launch_angle=-10.0, source_depth=100.0) + + +@pytest.fixture +def simple_rayfan(): + """Small RayFan of 3 synthetic Rays — no ray tracing needed.""" + rays = [ + _make_ray(launch_angle=-5.0, source_depth=100.0, n_bottom=0), + _make_ray(launch_angle=5.0, source_depth=150.0, n_bottom=1), + _make_ray(launch_angle=-10.0, source_depth=200.0, n_bottom=0), + ] + return RayFan(rays) diff --git a/tests/fixtures/munk_regression.npz b/tests/fixtures/munk_regression.npz new file mode 100644 index 0000000000000000000000000000000000000000..58b6fdbe07151247141997f31279501b72b84eac GIT binary patch literal 7586 zcmdU!c{o*F|NoEiAaOEfNRFvv9y7D|qHaSW5$=)kU zrcB8YA(ApgR7mmLp8M|p-p_MAzw7t>hCjZ)eO>!{oxRWgtiASnulHV`Ykw?uGcd}N z{@6rGo_X2|YkzLcBp#BFw}OY~IXV)BWWn;sHNO~dzTeEknvQgabb;pN?0wXWrcI^k zx@gkms5F;TUOrw99(Jd^oSfIs8#=gqI}_)IjWkHv>Y{n`hU8y{n=#_ z_nNO=-}iCUcvmb-^3qpybskB%c$Nh28nl0BpJf8ePAr^r4jX`M&rzmV$0*=SwNZ=V z7H+_5E6Ve!jTdx?1txa{@dE*J)3B(j5V-91AbjP$Feokqj@QkpVB1!Alc`&xz@BlJ zc_5282#e&py6Kn%7}P)e*0Mno_}@O5{Y8=nV%=iT6knDC;$qu5)yAbkcR){s>rPoP zJer&?mM;fBb=owU(J6q(mqmXYI-&>!qg+2t999AXzXyqPu`2`Dm%eZ2s+EEF3y$fl z$5en?8`XJ)M-^1$fxW9wRDqh`efLR^?cn2kj_-Gb)Bu@&InBIL4ODixZ}0O_2k+jU z@82h?0fg+bE1$p6092Q725?>z`1OLhr`xo^b#}XCyO&y^L+<1>e~>mv$Pz4EmeTX1l~1A{7lb?TqxLGN;bgq*n|7$j5X0`@C`zQ*X$ zb_QjDh~Ly^D^mue#<#D_*{gu`!GbIe4pq>ld{IdLu`2j+D>R(*_;#RTld)-Zvl@6G z7iE9*nHsn=e>Tv?TOB;U^t7Bx=+VtBdmv7+q~3L@?f(au@|Ft4`vBYx|@^tikAz}}{()&gzNKl-RjJRm%?Ie_r5u%jJS&06$2q3?r@+ z0x8`^?FpZR!Pcd>#PwA0=HhJZ*gaA3;pfUWC3bNzdyf3r-&F!c{@6ssOM=7%q4&I^ zG$3pxL9`X9l~aH_5Gcij%ctqqv-`P9Jg z4bLr_o~nVBQnU0~FLiK$+C9!o=y9Fvkdkh@2Jj6VJ?|Kx2_i2Xe|A+?3vhm%=e*ai z1wLh{jZo^|g zg=xTHT-vqZoD?`UK5+BID`~KH|K@RPWmz!RQl92@O%Bxgy$^l#Q6B7GJK>{Zs0ij* z7OLI%DgpHmrQctzDS?3!y7s+ym4PzHROEoI3fNG%Byw_tDqs~C(_E`o1@0Mb1>uAq zB8KQ31L2RPSNbw8gg>MPXQqn@e;f^Wb!#E?kP!cUXEUKktm^q~w+KDR&PBJp2t9`S z;$jpCJ&x>hbekmn5igMcyo~V2z09c7afCmB**BJ0!XFh46S`%DKccgoYrhiyD5nZw z2MB-kZOt6+B>bTuEv>fhj~^R1i?92`V(7@{zx(6!gD|ftV*Mz+TeVSk`77G+am3}F z9|=5oezIO^feF;cUy4yWu>r7tlU+OFNdaZjs{Iqv++gV=6Gdy37d&oi(s4-T2dR`d zD!zM!fFaX-U^L@api|$@SD0;f2bGm?>QB>X`)K>xkt3Gu#y z*Lj}^(f~h0!N@IdDM0%1x-q9!8c03p*;XMb3lw*A2!FpQ2Yl}>dk6Q)1LM_xS6dZD z@Zi83WR+8Ib%2jQ!z+hqT@bLPO7n}h z9taMZQ9n7M2WE%4GXxS4KzZebYU>m0NAX6_Z=Xyt*tQ)8bhH>+Jty&R*uukDL+2CjZjJTc2UnCp=gRAyv71dsZ>~DUG3-AWRqc`0L zQ&K|V)*XlUS3JK0t2}G0di^7yiKjbuCN3J5yhyeD;S~#W!r`%;JMpl8@mPP#$3#e# zusvHNk^+@_9@{S&rosFRukO%UXTZ_klj2`2vY=M`7TRUr9BA09aY$i42Nu;;WgV!^ zg?(@Cw@IAIgEu?kj!g^Y!vWftIoZ4Ua9ipj^+34-SeK?%mYP%mJGz64^`{FU?~MZi zrt*a_kgXkQxws_` zN=DpuByWm^XH*$C_j^ zFK9==hx%eU=9CE76?>pI)i?sSu}jSEeH#JYWjo7QUq!+^H`^E|p=@I%Sih(A>iEkp*V`0M=Nk8_TanQ6}&tK6t9_rPwm$RNrfI4+U)vAe!u+Am0y7zVx zGz(2GEw4(3ubZt{1glb@dEhm^p`uinmiYGl-jFmHY;}=e%{U!0C|yo0{+15689ir} zipqe`AK>a|$eFPCat=%Ou}m0*;pJ1eGoh@fNX*d7On7Q}@uuQ%Cj2())7tzZ6Os=x zP};LIAs~aW>-#gI)PwdLiL)8-4#Td=0*?%s?U(KYU#7z|TK^YN_Dj9+chK-li z33h(2Qr|(aRLmsCg5ai3@86FQ%*k-x;Z9IYVe6Lj1Q+()&k7@$oL3l?KyXQAsW^w= zUPpH9HbL(zFtLK*i%0bg7(u(12P01jh8vt=ZzibRVt48VK@C;9vvk`*f4KA801ihL18n+QV56Ibc7m-93&G3;)56lj=?O}mz3jX`&c%_%BJ146m{d(iu)puAEEBWB`(?n{ht3zA7BgTZeMNHcp-kv{3YF^5&xB2l;^E`%nXsSc{o4Uz z{^{8)w{mr5Ldl%5M#Y<%(3Xd}aoHggGOK^yJinF!^?CGfzYNKMdc`i=Q>N149KCnT zMV)jQemz}{$0rS9ua_2elT+dA2&)osBLxa#nvy+5$#6_7rc*IB3D&HB$|Q#*!u^fx zQ8%3vAk)R7(0Su{Xc%qaZ>1auTk@GF8YE)jo{-?46zLc!A5{`k7!?ga&)zFgelE%UBFC$>RNWan@{Rqe#MK861@V};~V%a+o0U1vFmw0tWz^dH2*2_5Y zeow@0zph2W(|#TnJHw;l{gQsMQd$hOz|@adh{Zxqwz;!Xig7TqwcM}VFdhnNs(x`k zngDsddIVYm6aUhqC?K)wTVxV!vmZU2f~CMOC*Ic`8A*ecuO}x?Ol82SZnCV}Ko)#v znIIr+l>-;(+5+Fn=fWC_h1tYHE-W*WRIt902Pw@co?(;^nfP>OHM{fSvqiPZ4$}h2 z>Rh@gTv`AV1}Y{eNrf=In{k$+TL||E8u0o$7Q(TwOeM@dg^=gWJ|oPt5V~`g?94n^ z2t#)r`B)-b2=|pfTw|Cd;&%oxV(|q~WPX@LmM(xx2e?pyTlw&bTIGoT=6sl*6f-pL znFq`A4nMm3Fc-2usx*H!l>>VlbFUP06;!8#{`R3_#I>6tct~F97oaE(xe!TS5tMQMU@%W5l_d90%cxlE$TevzuzO#CWDp0!_XX1A8wfM$|hkW=M z>$t#+pE;Dy|LqA6{?eb`uiTv*$Hr&PMpicAjTLrr-M$oDZv1}c1hx@BW2QU*(~k{r z(m5ylNQ8{fXNla*d&q#RFOVoEJ4txBDe{9Mb{X54AOA`3+B_Dkn|{`qegea05TQe` z6Px7w>YB~>46#=))-!h-Lu5Ddf3mEbLsBl>l8bI!M24QczmiqIiVP`_2t@bNqsh0^ zo!bnU(dkap0j(TXl!sjYHhPL3rBzH6lyFi|`NZAkrfHi{;kY#w$|i2~g0!s#n+Fei zP+VX+sF)W$eA^P7uI57-IlK-&a9}~Nu&&uhZeu}@iA_dUpCqHE+8@s!muE())pV0K zc8n<1TXX9hMKdy8@M28A7F`HTe)+U8&PU&2Ckt?xK*UBiSYlt)Az>G37^>g4QOjJSYSl^3rQ zGwv*`@xd*IjPsnhWE8B6;Q7zEuLTLR;D^eFOnFX|@dLv1n`S(jaZT>7 z`4vef9K(}v({u(rkHnn+>NOoMtza0+FtCO#H^=uJsa(OVq7AjH_Ag={?~F|^KAXqz zj@~1UqSM&%{nj+0OQV<|!pw9(yc6pfPR)Lnn}d{^h}WOI`wGFBIqgPb-XT}tm#wf% zenPgw3S%?4fb5~u^*Cj;gvjA}K3*9=5dXrdR7-0T+7v%!GOS0B@_d;JdAPuUK6f^_5e z)LseW!LA&4h*LW#f$9G6bxh(=z{XZSk}hSaVlJWO!O1-um^Mi*kCvx{z2PxUt(-wH z!&HW$`EGzk_n{i^-|1uTU;Mb$k!yglI~LX0v+Tfv%db7H!zgAUS)S1hb+JidvzE9mT9_2gCT)^I4ZFw9^{g{o2~(?b z>ZDwuVb-fG<2U()FltM;zQeWewf=M`{Kf}2AUZpG>>i7VB7XBej-Qugkat@@N8XoF zK^|RLatyhoj&K!4de-`CBU4^)NG8*I2ygJ0MIm>9WZvyK^&wFoITmcH0F4 zQ8Q;?YoWPHo5T3S&WGcJ0ZzXP{s z!jE3K`8>6k0iTW;tFqWikB@dHi+oZi;diu7eqGr21M53)6Sd|0BG$9@*;%sZ0;WB@ zQpgZAgDq5>XcT`K!*Y7qh16}ju{iZ;VaM(;%$Ps_js8#h^pZ zd75M&tjM>3TxDH~xSYR;jBeHq3m_5mZ|4_^F@=O4v!OhB^_>pg#dg1b<{<;xT57RB z1z|#aWjHp!xWtV9$nIo0A4o>=m%WFh#8}Yt3eDy*sf@o|zl(zm_TTD0&la zKl%A{5(5Q~6{2>=b+O|cpIs}>y26TI^Brq?D8-D^yhC{fYw7WWdb&4cqgOHiE!jTz z;}$V$EAM8xs5wl8pKH&JAETIi37YPGv>wZVl@ZgaT?ngOe&vYA1k!wle9F3h9%=R~ z59w%FMtaG5AxlR{=$?)pT;whWl$Ki8^lLsw65or_~bU$7$^YI--K3pbksI)W)E zS3rcnDVYn68h-ZC`7Ado`>`jZvV{k=9&6N|Ug1R_=#?t$_{oR%DrFa~KG}@!+O;FL zU6UVG8*#yMO!?8TT`Uq+?>3|Su`>O_@A=SAT2c4*sPm!)Ls7H(kP9_*DiYTI#(~Do z3})RPWkxg0_uOgbTR}o{OPm%}29Tjkn`vY1Q)q;R#Vh%D<}a5~ zaHVwS_f7Y?aayQ+;G+><{Gam>7sge&UwOy|qAufM=XlEJzlt+A|9rds;=eD>xc#Ta zng6TG3{lbft2py#xBjR2Up=I^ub0bT8P~6e`^Eo1AM#|kpT;jG+14v;uFgIV|C`VF zA6MD_9|=i?e<(=#lfU$&fIogwiLZBb^zsbjM7MtM$He{DnDuXd z|32m_(aIlR0RJ^+JwN~VF_cXt(jVFSzs9V4X?+f@?;F;;^XC>smm)y=yR-iDzyQJZ tz`thezXq;r{_Bo^_eKz%qYyC~f6V}k-Aw-&EM{WN42h94+_K)a{{gq;aby4h literal 0 HcmV?d00001 diff --git a/tests/test_environment.py b/tests/test_environment.py new file mode 100644 index 0000000..a639236 --- /dev/null +++ b/tests/test_environment.py @@ -0,0 +1,180 @@ +""" +Tests for pygenray.environment: munk_ssp, OceanEnvironment2D, eflat, eflatinv. +""" +import numpy as np +import pytest +import xarray as xr +from matplotlib import pyplot as plt + +from pygenray.environment import ( + OceanEnvironment2D, + eflat, + eflatinv, + munk_ssp, +) + + +# --------------------------------------------------------------------------- +# munk_ssp +# --------------------------------------------------------------------------- + +class TestMunkSSP: + def test_output_shape_matches_input(self): + z = np.arange(0, 5000, 10) + c = munk_ssp(z) + assert c.shape == z.shape + + def test_minimum_at_sofar_depth(self): + sofar = 1300.0 + z = np.arange(0, 6000, 1) + c = munk_ssp(z, sofar_depth=sofar) + # Minimum should be at the SOFAR depth + assert z[np.argmin(c)] == pytest.approx(sofar, abs=2.0) + + def test_default_params_near_1500_at_sofar(self): + sofar = 1300.0 + c_sofar = munk_ssp(np.array([sofar])) + assert c_sofar[0] == pytest.approx(1500.0, abs=5.0) + + def test_scalar_input(self): + c = munk_ssp(np.array([0.0])) + assert c.shape == (1,) + + +# --------------------------------------------------------------------------- +# OceanEnvironment2D construction +# --------------------------------------------------------------------------- + +class TestOceanEnvironment2DConstruction: + def test_default_init_attributes_exist(self): + env = OceanEnvironment2D() + for attr in ('sound_speed', 'bathymetry', 'dcdz', 'bottom_angle', + 'bottom_angle_interp'): + assert hasattr(env, attr), f"Missing attribute: {attr}" + + def test_default_sound_speed_is_2d(self): + env = OceanEnvironment2D() + assert env.sound_speed.ndim == 2 + assert set(env.sound_speed.dims) == {'range', 'depth'} + + def test_default_flat_earth_attributes_exist(self): + env = OceanEnvironment2D(flat_earth_transform=True) + assert hasattr(env, 'sound_speed_fe') + assert hasattr(env, 'bathymetry_fe') + + def test_flat_earth_false_no_fe_attributes(self): + env = OceanEnvironment2D(flat_earth_transform=False) + assert not hasattr(env, 'sound_speed_fe') + assert not hasattr(env, 'bathymetry_fe') + + def test_custom_1d_sound_speed(self): + z = np.arange(0.0, 3000.0, 10.0) + c_vals = munk_ssp(z) + ssp = xr.DataArray(c_vals, dims=['depth'], coords={'depth': z}) + bathy = xr.DataArray( + np.ones(20) * 4000.0, dims=['range'], + coords={'range': np.linspace(0, 50e3, 20)} + ) + env = OceanEnvironment2D(sound_speed=ssp, bathymetry=bathy, + flat_earth_transform=False) + assert env.sound_speed.ndim == 1 + assert 'depth' in env.sound_speed.dims + + def test_custom_2d_sound_speed(self): + z = np.arange(0.0, 3000.0, 50.0) + r = np.linspace(0.0, 50e3, 20) + c_2d = np.outer(np.ones(len(r)), munk_ssp(z)) + ssp = xr.DataArray(c_2d, dims=['range', 'depth'], + coords={'range': r, 'depth': z}) + env = OceanEnvironment2D(sound_speed=ssp, flat_earth_transform=False) + assert env.sound_speed.ndim == 2 + + def test_custom_bathymetry_stored(self): + bathy_vals = np.ones(20) * 3500.0 + r = np.linspace(0.0, 50e3, 20) + bathy = xr.DataArray(bathy_vals, dims=['range'], coords={'range': r}) + env = OceanEnvironment2D(bathymetry=bathy, flat_earth_transform=False) + np.testing.assert_array_equal(env.bathymetry.values, bathy_vals) + + # --- invalid inputs --- + + def test_sound_speed_not_dataarray_raises_type_error(self): + with pytest.raises(TypeError): + OceanEnvironment2D(sound_speed=np.ones(100)) + + def test_sound_speed_3d_raises_value_error(self): + da = xr.DataArray( + np.ones((5, 10, 20)), + dims=['range', 'depth', 'extra'], + coords={'range': np.arange(5), 'depth': np.arange(10), + 'extra': np.arange(20)} + ) + with pytest.raises(ValueError): + OceanEnvironment2D(sound_speed=da) + + def test_sound_speed_missing_depth_dim_raises_value_error(self): + da = xr.DataArray(np.ones(50), dims=['range'], + coords={'range': np.arange(50)}) + with pytest.raises(ValueError): + OceanEnvironment2D(sound_speed=da) + + def test_2d_sound_speed_missing_range_dim_raises_value_error(self): + da = xr.DataArray( + np.ones((10, 20)), + dims=['depth', 'extra'], + coords={'depth': np.arange(10), 'extra': np.arange(20)} + ) + with pytest.raises(ValueError): + OceanEnvironment2D(sound_speed=da) + + def test_bathymetry_not_dataarray_raises_type_error(self): + with pytest.raises(TypeError): + OceanEnvironment2D(bathymetry=np.ones(50)) + + def test_bathymetry_missing_range_dim_raises_value_error(self): + da = xr.DataArray(np.ones(50), dims=['depth'], + coords={'depth': np.arange(50)}) + with pytest.raises(ValueError): + OceanEnvironment2D(bathymetry=da) + + +# --------------------------------------------------------------------------- +# eflat / eflatinv round-trip +# --------------------------------------------------------------------------- + +class TestEflat: + LAT = 35.0 + + def test_depth_roundtrip(self): + dep = np.array([100.0, 500.0, 1000.0, 2000.0, 4000.0]) + depf, _ = eflat(dep, self.LAT) + dep_rec, _ = eflatinv(depf, np.array([self.LAT])) + np.testing.assert_allclose(dep_rec, dep, atol=1.0, + err_msg="Depth round-trip outside 1 m tolerance") + + def test_sound_speed_roundtrip(self): + dep = np.array([100.0, 500.0, 1000.0, 2000.0]) + cs = np.array([1500.0, 1490.0, 1480.0, 1510.0]) + depf, csf = eflat(dep, self.LAT, cs) + _, cs_rec = eflatinv(depf, np.array([self.LAT]), csf) + np.testing.assert_allclose(cs_rec, cs, rtol=1e-4, + err_msg="Sound speed round-trip outside 0.01% tolerance") + + def test_eflat_increases_depth(self): + """Flat-earth transformation should increase effective depths.""" + dep = np.array([100.0, 1000.0, 3000.0]) + depf, _ = eflat(dep, self.LAT) + assert np.all(depf > dep) + + +# --------------------------------------------------------------------------- +# OceanEnvironment2D.plot smoke test +# --------------------------------------------------------------------------- + +class TestOceanEnvironment2DPlot: + def test_plot_runs_without_error(self): + env = OceanEnvironment2D() + fig, ax = plt.subplots() + plt.sca(ax) + env.plot() + plt.close('all') diff --git a/tests/test_physics.py b/tests/test_physics.py new file mode 100644 index 0000000..499314f --- /dev/null +++ b/tests/test_physics.py @@ -0,0 +1,324 @@ +""" +Physics correctness tests for pygenray ray tracing. + +All tests use flatearth=False to avoid flat-earth correction complicating +the comparison with analytical solutions. +""" +import os +import pathlib + +import numpy as np +import pytest +import xarray as xr + +from pygenray.environment import OceanEnvironment2D, munk_ssp +from pygenray.launch_rays import shoot_ray, shoot_rays + +FIXTURE_DIR = pathlib.Path(__file__).parent / 'fixtures' + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _const_c_env(c0=1500.0, z_max=5000.0, r_max=100e3, + bathy_depth=4500.0, nz=200, nr=20): + """Build a range-independent constant sound-speed environment.""" + z = np.linspace(0.0, z_max, nz) + r = np.linspace(0.0, r_max, nr) + c_2d = np.full((nr, nz), c0) + ssp = xr.DataArray(c_2d, dims=['range', 'depth'], + coords={'range': r, 'depth': z}) + bathy = xr.DataArray(np.full(nr, bathy_depth), dims=['range'], + coords={'range': r}) + return OceanEnvironment2D(sound_speed=ssp, bathymetry=bathy, + flat_earth_transform=False) + + +def _linear_gradient_env(c0=1500.0, g=0.05, z_max=5000.0, r_max=100e3, + bathy_depth=4500.0, nz=500, nr=50): + """Build a range-independent linear-gradient environment: c(z) = c0 + g*z.""" + z = np.linspace(0.0, z_max, nz) + r = np.linspace(0.0, r_max, nr) + c_1d = c0 + g * z + c_2d = np.outer(np.ones(nr), c_1d) + ssp = xr.DataArray(c_2d, dims=['range', 'depth'], + coords={'range': r, 'depth': z}) + bathy = xr.DataArray(np.full(nr, bathy_depth), dims=['range'], + coords={'range': r}) + return OceanEnvironment2D(sound_speed=ssp, bathymetry=bathy, + flat_earth_transform=False) + + +def _munk_env(r_max=100e3, nr=50, nz=600, bathy_depth=5000.0): + """Build a range-independent Munk-profile environment.""" + z = np.linspace(0.0, 6000.0, nz) + r = np.linspace(0.0, r_max, nr) + c_1d = munk_ssp(z) + c_2d = np.outer(np.ones(nr), c_1d) + ssp = xr.DataArray(c_2d, dims=['range', 'depth'], + coords={'range': r, 'depth': z}) + bathy = xr.DataArray(np.full(nr, bathy_depth), dims=['range'], + coords={'range': r}) + return OceanEnvironment2D(sound_speed=ssp, bathymetry=bathy, + flat_earth_transform=False) + + +# --------------------------------------------------------------------------- +# A. Snell's invariant in range-independent constant-c medium +# --------------------------------------------------------------------------- + +class TestSnellInvariant: + """ + In a constant sound-speed medium (dc/dz = 0), the ray-parameter + p = sin(θ)/c satisfies dp/dx = 0 exactly. The stored value ray.p = -p_ode + should therefore be constant to within ODE integration error. + """ + + @pytest.mark.parametrize("user_angle", [-5.0, -10.0, -15.0]) + def test_p_constant_along_ray(self, user_angle): + c0 = 1500.0 + env = _const_c_env(c0=c0) + ray = shoot_ray(200.0, 0.0, user_angle, 30e3, 60, + env, rtol=1e-9, flatearth=False, debug=False) + assert ray is not None, "shoot_ray returned None unexpectedly" + # |p| should be constant; sign may flip at boundary reflections + abs_p = np.abs(ray.p) + p_std = np.std(abs_p) + p_mean = np.mean(abs_p) + assert p_std / p_mean < 1e-5, ( + f"|p| not constant in constant-c medium: std/mean = {p_std/p_mean:.2e}" + ) + + +# --------------------------------------------------------------------------- +# B. Constant sound speed — straight-line rays +# --------------------------------------------------------------------------- + +class TestConstantSSPStraightLine: + """ + In a homogeneous medium the ray paths are straight lines. + Analytical travel time: t = R / (c · cos θ) + Analytical final depth: z = z0 + R · tan θ (positive downward, ODE convention) + Stored convention: ray.z = -(ODE z) + """ + + def test_travel_time_analytical(self): + c0 = 1500.0 + user_angle = -10.0 # downward + theta_ode = abs(user_angle) # internal ODE angle (degrees) + z0 = 200.0 # source depth [m] + R = 20e3 # receiver range [m] + env = _const_c_env(c0=c0, r_max=R + 1e3) + + ray = shoot_ray(z0, 0.0, user_angle, R, 50, + env, rtol=1e-9, flatearth=False, debug=False) + assert ray is not None + + t_analytical = R / (c0 * np.cos(np.radians(theta_ode))) + assert abs(ray.t[-1] - t_analytical) / t_analytical < 1e-3 + + def test_final_depth_analytical(self): + c0 = 1500.0 + user_angle = -10.0 + theta_ode = abs(user_angle) + z0 = 200.0 + R = 20e3 + env = _const_c_env(c0=c0, r_max=R + 1e3) + + ray = shoot_ray(z0, 0.0, user_angle, R, 50, + env, rtol=1e-9, flatearth=False, debug=False) + assert ray is not None + + z_ode_end = z0 + R * np.tan(np.radians(theta_ode)) + z_expected = -z_ode_end # stored sign convention + assert abs(ray.z[-1] - z_expected) / abs(z_expected) < 1e-3 + + def test_p_constant_in_const_c(self): + c0 = 1500.0 + user_angle = -10.0 + theta_ode = abs(user_angle) + z0 = 200.0 + R = 20e3 + env = _const_c_env(c0=c0, r_max=R + 1e3) + + ray = shoot_ray(z0, 0.0, user_angle, R, 50, + env, rtol=1e-9, flatearth=False, debug=False) + assert ray is not None + + p_expected = -np.sin(np.radians(theta_ode)) / c0 # stored sign + np.testing.assert_allclose(ray.p, p_expected, + rtol=1e-5, atol=0, + err_msg="p not constant in homogeneous medium") + + +# --------------------------------------------------------------------------- +# C. Linear gradient — turning depth +# --------------------------------------------------------------------------- + +class TestLinearGradientTurningDepth: + """ + For c(z) = c0 + g·z the Hamiltonian is conserved: + H = sqrt(1/c(z)² - p_ode²) = sqrt(1/c_source² - p0²) + At the turning point (ray horizontal, dz/dx = 0) we have p_ode = 0, so: + c(z_turn) = 1 / sqrt(1/c_source² - p0²) + where p0 = sin(θ_ode) / c_source. + + Substituting: c_source = c0 + g·z_source, p0 = sin(θ)/c_source + c(z_turn) = c_source / cos(θ) + z_turn = (c_source/cos(θ) − c0) / g + """ + + C0 = 1500.0 + G = 0.05 # m/s/m + Z_SRC = 200.0 # source depth [m] + THETA = 20.0 # degrees downward (user passes -20) + + def _z_turn_analytical(self): + c_source = self.C0 + self.G * self.Z_SRC + return (c_source / np.cos(np.radians(self.THETA)) - self.C0) / self.G + + def test_turning_depth_approx(self): + z_turn = self._z_turn_analytical() + env = _linear_gradient_env(c0=self.C0, g=self.G) + + ray = shoot_ray(self.Z_SRC, 0.0, -self.THETA, 80e3, 400, + env, rtol=1e-9, flatearth=False, debug=False) + assert ray is not None + + # ray.z uses sign convention: negative = deep + z_turn_numerical = -np.min(ray.z) # max depth reached + assert abs(z_turn_numerical - z_turn) < 50.0, ( + f"Turning depth: expected {z_turn:.1f} m, got {z_turn_numerical:.1f} m" + ) + + def test_hamiltonian_conserved_linear_gradient(self): + """The Hamiltonian H = sqrt(1/c(z)² − p_ode²) is conserved.""" + env = _linear_gradient_env(c0=self.C0, g=self.G) + + ray = shoot_ray(self.Z_SRC, 0.0, -self.THETA, 80e3, 400, + env, rtol=1e-9, flatearth=False, debug=False) + assert ray is not None + + z_ode = -ray.z # positive depth (ODE convention) + p_ode = -ray.p # positive for downward ray + c_along = self.C0 + self.G * z_ode + + H = np.sqrt(1.0 / c_along**2 - p_ode**2) + H_std = np.std(H) + H_mean = np.mean(H) + assert H_std / H_mean < 1e-4, ( + f"Hamiltonian not conserved: std/mean = {H_std/H_mean:.2e}" + ) + + +# --------------------------------------------------------------------------- +# D. Hamiltonian conservation in range-independent Munk profile +# --------------------------------------------------------------------------- + +class TestMunkHamiltonianConservation: + """ + In any range-independent environment the Hamiltonian + H = sqrt(1/c(z)² − p_ode²) + is an exact invariant of the ray equations. Numerical integration + should preserve it to within ODE tolerance. + """ + + @pytest.mark.parametrize("user_angle", [-5.0, -10.0, -15.0]) + def test_hamiltonian_conserved_munk(self, user_angle): + env = _munk_env(r_max=100e3) + z_src = 1000.0 # below SOFAR axis — ray stays in the water column + + ray = shoot_ray(z_src, 0.0, user_angle, 100e3, 200, + env, rtol=1e-9, flatearth=False, debug=False) + assert ray is not None, "shoot_ray returned None; ray may have exited domain" + + z_ode = -ray.z + p_ode = -ray.p + c_along = munk_ssp(z_ode) + + arg = 1.0 / c_along**2 - p_ode**2 + # Clamp small negatives from floating-point near turning point + arg = np.clip(arg, 0.0, None) + H = np.sqrt(arg) + + # Exclude points at/near the turning point where H → 0 + mask = H > 1e-6 / 1500.0 + if mask.sum() < 5: + pytest.skip("Too few valid points away from turning point") + + H_valid = H[mask] + H_std = np.std(H_valid) + H_mean = np.mean(H_valid) + assert H_std / H_mean < 1e-3, ( + f"Hamiltonian not conserved in Munk profile: std/mean = {H_std/H_mean:.2e}" + ) + + +# --------------------------------------------------------------------------- +# E. Regression / golden-file tests +# --------------------------------------------------------------------------- + +def _regenerate_fixture(): + """Run shoot_rays and save regression fixture. Call manually to regenerate.""" + env = _munk_env(r_max=50e3, nr=30, nz=400) + angles = [-8.0, -4.0, 0.0, 4.0, 8.0] + rf = shoot_rays(1300.0, 0.0, angles, 50e3, 50, env, + n_processes=1, debug=False, flatearth=False) + FIXTURE_DIR.mkdir(parents=True, exist_ok=True) + np.savez( + FIXTURE_DIR / 'munk_regression.npz', + ts=rf.ts, zs=rf.zs, ps=rf.ps, + n_botts=rf.n_botts, n_surfs=rf.n_surfs, + thetas=rf.thetas, + ) + return rf + + +class TestMunkRegression: + """ + Golden-file regression: compare shoot_rays output against a saved fixture. + + To regenerate the fixture (e.g., after intentional physics changes): + python -c "from tests.test_physics import _regenerate_fixture; _regenerate_fixture()" + or run pytest with --regenerate-physics flag. + """ + + FIXTURE = FIXTURE_DIR / 'munk_regression.npz' + ANGLES = [-8.0, -4.0, 0.0, 4.0, 8.0] + + def _run_rays(self): + env = _munk_env(r_max=50e3, nr=30, nz=400) + return shoot_rays(1300.0, 0.0, self.ANGLES, 50e3, 50, env, + n_processes=1, debug=False, flatearth=False) + + def test_regression(self, request): + regenerate = request.config.getoption('--regenerate-physics', + default=False) + if regenerate or not self.FIXTURE.exists(): + rf = self._run_rays() + self.FIXTURE.parent.mkdir(parents=True, exist_ok=True) + np.savez( + self.FIXTURE, + ts=rf.ts, zs=rf.zs, ps=rf.ps, + n_botts=rf.n_botts, n_surfs=rf.n_surfs, + thetas=rf.thetas, + ) + if regenerate: + pytest.skip("Fixture regenerated; skipping comparison") + # First-run: fixture just created, nothing to compare + return + + ref = np.load(self.FIXTURE) + rf = self._run_rays() + + np.testing.assert_allclose(rf.ts, ref['ts'], atol=1e-6, + err_msg="Travel times changed") + np.testing.assert_allclose(rf.zs, ref['zs'], atol=1e-6, + err_msg="Depths changed") + np.testing.assert_allclose(rf.ps, ref['ps'], atol=1e-6, + err_msg="Ray parameters changed") + np.testing.assert_array_equal(rf.n_botts, ref['n_botts']) + np.testing.assert_array_equal(rf.n_surfs, ref['n_surfs']) + + diff --git a/tests/test_ray_objects.py b/tests/test_ray_objects.py new file mode 100644 index 0000000..509dd68 --- /dev/null +++ b/tests/test_ray_objects.py @@ -0,0 +1,278 @@ +""" +Tests for pygenray.ray_objects: Ray, RayFan. +""" +import numpy as np +import pytest +import scipy.io +from matplotlib import pyplot as plt + +from pygenray.ray_objects import Ray, RayFan + + +# --------------------------------------------------------------------------- +# Ray +# --------------------------------------------------------------------------- + +class TestRay: + N = 10 + R = 10000.0 + + def _make_ray(self, launch_angle=-10.0, source_depth=100.0, + n_bottom=0, n_surface=0): + r = np.linspace(0.0, self.R, self.N) + t = r / 1500.0 + z_ode = np.linspace(source_depth, source_depth + self.R * 0.01, self.N) + p_ode = np.ones(self.N) * np.sin(np.radians(abs(launch_angle) + 1e-3)) / 1500.0 + y = np.vstack([t, z_ode, p_ode]) + return Ray(r=r, y=y, n_bottom=n_bottom, n_surface=n_surface, + launch_angle=launch_angle, source_depth=source_depth), y + + def test_attribute_shapes(self): + ray, _ = self._make_ray() + assert ray.r.shape == (self.N,) + assert ray.t.shape == (self.N,) + assert ray.z.shape == (self.N,) + assert ray.p.shape == (self.N,) + + def test_z_sign_convention(self): + """ray.z should equal -y[1,:].""" + ray, y = self._make_ray() + np.testing.assert_array_equal(ray.z, -y[1, :]) + + def test_p_sign_convention(self): + """ray.p should equal -y[2,:].""" + ray, y = self._make_ray() + np.testing.assert_array_equal(ray.p, -y[2, :]) + + def test_launch_angle_stored(self): + ray, _ = self._make_ray(launch_angle=-15.0) + assert ray.launch_angle == pytest.approx(-15.0) + + def test_source_depth_stored(self): + ray, _ = self._make_ray(source_depth=250.0) + assert ray.source_depth == pytest.approx(250.0) + + def test_optional_launch_angle_not_set(self): + r = np.linspace(0.0, self.R, self.N) + t = r / 1500.0 + y = np.vstack([t, np.ones(self.N) * 100.0, np.ones(self.N) * 0.1]) + ray = Ray(r=r, y=y, n_bottom=0, n_surface=0) + assert not hasattr(ray, 'launch_angle') + + def test_optional_source_depth_not_set(self): + r = np.linspace(0.0, self.R, self.N) + t = r / 1500.0 + y = np.vstack([t, np.ones(self.N) * 100.0, np.ones(self.N) * 0.1]) + ray = Ray(r=r, y=y, n_bottom=0, n_surface=0) + assert not hasattr(ray, 'source_depth') + + def test_n_bottom_n_surface_stored(self): + ray, _ = self._make_ray(n_bottom=3, n_surface=1) + assert ray.n_bottom == 3 + assert ray.n_surface == 1 + + def test_plot_smoke(self): + ray, _ = self._make_ray() + plt.figure() + ray.plot() + plt.close('all') + + +# --------------------------------------------------------------------------- +# RayFan +# --------------------------------------------------------------------------- + +class TestRayFan: + M = 3 + N = 10 + R = 10000.0 + + def _make_rays(self, M=None, N=None, R=None): + M = M or self.M + N = N or self.N + R = R or self.R + rays = [] + for i in range(M): + r = np.linspace(0.0, R, N) + theta = float(-5 + i * 5) + t = r / 1500.0 + z_ode = np.linspace(100.0 + i * 50, 200.0 + i * 50, N) + p_ode = np.ones(N) * np.sin(np.radians(abs(theta) + 1e-3)) / 1500.0 + y = np.vstack([t, z_ode, p_ode]) + rays.append( + Ray(r=r, y=y, n_bottom=i % 2, n_surface=0, + launch_angle=theta, source_depth=100.0 + i * 50) + ) + return rays + + # --- Construction --- + + def test_shapes(self, simple_rayfan): + rf = simple_rayfan + assert rf.thetas.shape == (self.M,) + assert rf.rs.shape == (self.M, self.N) + assert rf.ts.shape == (self.M, self.N) + assert rf.zs.shape == (self.M, self.N) + assert rf.ps.shape == (self.M, self.N) + assert rf.n_botts.shape == (self.M,) + assert rf.n_surfs.shape == (self.M,) + assert rf.source_depths.shape == (self.M,) + + def test_ray_ids_set_on_construction(self, simple_rayfan): + assert hasattr(simple_rayfan, 'ray_ids') + assert len(simple_rayfan.ray_ids) == self.M + + # --- compute_rayids --- + + def test_compute_rayids_returns_strings(self, simple_rayfan): + simple_rayfan.compute_rayids() + assert all(isinstance(rid, str) for rid in simple_rayfan.ray_ids) + + def test_compute_rayids_length(self, simple_rayfan): + simple_rayfan.compute_rayids() + assert len(simple_rayfan.ray_ids) == len(simple_rayfan.thetas) + + # --- __len__ --- + + def test_len(self, simple_rayfan): + assert len(simple_rayfan) == self.M + + # --- __getitem__ int --- + + def test_getitem_int_returns_ray(self, simple_rayfan): + ray = simple_rayfan[0] + assert isinstance(ray, Ray) + + def test_getitem_int_correct_index(self, simple_rayfan): + ray = simple_rayfan[1] + np.testing.assert_array_equal(ray.r, simple_rayfan.rs[1]) + + def test_getitem_negative_int(self, simple_rayfan): + ray = simple_rayfan[-1] + assert isinstance(ray, Ray) + np.testing.assert_array_equal(ray.r, simple_rayfan.rs[-1]) + + def test_getitem_out_of_bounds_raises_index_error(self, simple_rayfan): + with pytest.raises(IndexError): + _ = simple_rayfan[100] + + # --- __getitem__ slice --- + + def test_getitem_slice_returns_rayfan(self, simple_rayfan): + result = simple_rayfan[0:2] + assert isinstance(result, RayFan) + assert len(result) == 2 + + def test_getitem_slice_correct_thetas(self, simple_rayfan): + result = simple_rayfan[1:] + np.testing.assert_array_equal(result.thetas, simple_rayfan.thetas[1:]) + + # --- __getitem__ bool mask --- + + def test_getitem_bool_mask_returns_rayfan(self, simple_rayfan): + mask = np.array([True, False, True]) + result = simple_rayfan[mask] + assert isinstance(result, RayFan) + assert len(result) == 2 + + def test_getitem_bool_mask_correct_subset(self, simple_rayfan): + mask = np.array([False, True, False]) + result = simple_rayfan[mask] + np.testing.assert_array_equal(result.thetas, simple_rayfan.thetas[1:2]) + + # --- __getitem__ int array --- + + def test_getitem_int_array_returns_rayfan(self, simple_rayfan): + idx = np.array([0, 2]) + result = simple_rayfan[idx] + assert isinstance(result, RayFan) + assert len(result) == 2 + np.testing.assert_array_equal(result.thetas, + simple_rayfan.thetas[np.array([0, 2])]) + + # --- __add__ --- + + def test_add_correct_length(self): + rays_a = self._make_rays(M=2) + rays_b = self._make_rays(M=3) + rf_a = RayFan(rays_a) + rf_b = RayFan(rays_b) + result = rf_a + rf_b + assert len(result) == 5 + + def test_add_rs_preserved(self): + rays_a = self._make_rays(M=2) + rays_b = self._make_rays(M=1) + rf_a = RayFan(rays_a) + rf_b = RayFan(rays_b) + result = rf_a + rf_b + # All rays should have the same range array + for i in range(len(result)): + np.testing.assert_array_equal(result.rs[i], rf_a.rs[0]) + + def test_add_incompatible_ranges_raises_value_error(self): + rays_a = self._make_rays(M=1, R=10000.0) + rays_b = self._make_rays(M=1, R=20000.0) + rf_a = RayFan(rays_a) + rf_b = RayFan(rays_b) + with pytest.raises(ValueError): + _ = rf_a + rf_b + + def test_add_non_rayfan_raises_type_error(self, simple_rayfan): + with pytest.raises(TypeError): + _ = simple_rayfan + 42 + + # --- save_mat --- + + def test_save_mat_creates_file(self, simple_rayfan, tmp_path): + path = str(tmp_path / 'test_rayfan.mat') + simple_rayfan.save_mat(path) + assert (tmp_path / 'test_rayfan.mat').exists() + + def test_save_mat_loadable(self, simple_rayfan, tmp_path): + path = str(tmp_path / 'test_rayfan.mat') + simple_rayfan.save_mat(path) + data = scipy.io.loadmat(path) + assert 'rayfan' in data + + def test_save_mat_contains_required_keys(self, simple_rayfan, tmp_path): + path = str(tmp_path / 'test_rayfan.mat') + simple_rayfan.save_mat(path) + data = scipy.io.loadmat(path) + rayfan = data['rayfan'] + # MATLAB struct stored as structured array; check dtype field names + expected_keys = {'thetas', 'xs', 'ts', 'zs', 'ps', + 'n_botts', 'n_surfs', 'source_depths'} + actual_keys = set(rayfan.dtype.names) + assert expected_keys <= actual_keys + + def test_save_mat_values_match(self, simple_rayfan, tmp_path): + path = str(tmp_path / 'test_rayfan.mat') + simple_rayfan.save_mat(path) + data = scipy.io.loadmat(path) + rayfan = data['rayfan'] + saved_thetas = rayfan['thetas'][0, 0].flatten() + np.testing.assert_allclose(saved_thetas, simple_rayfan.thetas, + atol=1e-10) + + # --- Plotting smoke tests --- + + def test_plot_ray_fan_smoke(self, simple_rayfan): + plt.figure() + simple_rayfan.plot_ray_fan() + plt.close('all') + + def test_plot_time_front_smoke(self, simple_rayfan): + plt.figure() + simple_rayfan.plot_time_front() + plt.close('all') + + def test_plot_time_front_include_lines_smoke(self, simple_rayfan): + plt.figure() + simple_rayfan.plot_time_front(include_lines=True) + plt.close('all') + + def test_plot_depth_v_angle_smoke(self, simple_rayfan): + plt.figure() + simple_rayfan.plot_depth_v_angle() + plt.close('all') From 569c6c6d0702f31fe35febc932f1a23eaad0fe40 Mon Sep 17 00:00:00 2001 From: John Ragland Date: Wed, 18 Mar 2026 15:47:00 -0400 Subject: [PATCH 2/2] reduced tolerence on munk --- tests/test_physics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_physics.py b/tests/test_physics.py index 499314f..62cae0d 100644 --- a/tests/test_physics.py +++ b/tests/test_physics.py @@ -314,9 +314,9 @@ def test_regression(self, request): np.testing.assert_allclose(rf.ts, ref['ts'], atol=1e-6, err_msg="Travel times changed") - np.testing.assert_allclose(rf.zs, ref['zs'], atol=1e-6, + np.testing.assert_allclose(rf.zs, ref['zs'], atol=0.1, err_msg="Depths changed") - np.testing.assert_allclose(rf.ps, ref['ps'], atol=1e-6, + np.testing.assert_allclose(rf.ps, ref['ps'], atol=0.1, err_msg="Ray parameters changed") np.testing.assert_array_equal(rf.n_botts, ref['n_botts']) np.testing.assert_array_equal(rf.n_surfs, ref['n_surfs'])