From 1f0f9f1bfc425cfa01e80befba765abcf6c55443 Mon Sep 17 00:00:00 2001 From: Philip Crotwell Date: Fri, 6 Sep 2024 17:21:39 -0400 Subject: [PATCH] allow multiple depths in single run, read depths from quakeml file arg if given. Remove duplicate Locations from stationxml. JSON has sourceDepths as list and each Arrival includes the source depth calculated from. Other refactor of JSON output code. --- docs/manual/taup_curve.html | 11 +- docs/manual/taup_find.html | 6 +- docs/manual/taup_path.html | 28 +-- docs/manual/taup_phase.html | 9 +- docs/manual/taup_pierce.html | 26 +-- docs/manual/taup_refltrans.html | 7 +- docs/manual/taup_setms3.html | 24 +-- docs/manual/taup_setsac.html | 6 +- docs/manual/taup_spikes.html | 16 +- docs/manual/taup_table.html | 11 +- docs/manual/taup_time.html | 14 +- docs/manual/taup_wavefront.html | 11 +- .../sphinx/source/cmdLineHelp/VERSION.json | 10 +- .../source/cmdLineHelp/taup_curve.usage | 11 +- .../sphinx/source/cmdLineHelp/taup_find.usage | 6 +- .../sphinx/source/cmdLineHelp/taup_path.usage | 28 +-- .../source/cmdLineHelp/taup_phase.usage | 9 +- .../source/cmdLineHelp/taup_pierce.usage | 26 +-- .../source/cmdLineHelp/taup_refltrans.usage | 7 +- .../source/cmdLineHelp/taup_setms3.usage | 24 +-- .../source/cmdLineHelp/taup_setsac.usage | 6 +- .../source/cmdLineHelp/taup_spikes.usage | 16 +- .../source/cmdLineHelp/taup_table.usage | 11 +- .../sphinx/source/cmdLineHelp/taup_time.usage | 14 +- .../source/cmdLineHelp/taup_wavefront.usage | 11 +- src/main/java/edu/sc/seis/TauP/Arrival.java | 3 + .../edu/sc/seis/TauP/RayCalculateable.java | 14 ++ src/main/java/edu/sc/seis/TauP/Scatterer.java | 5 + .../java/edu/sc/seis/TauP/XYPlotOutput.java | 3 +- .../TauP/cmdline/TauP_AbstractPhaseTool.java | 77 ++++++-- .../TauP/cmdline/TauP_AbstractRayTool.java | 129 ++++++++----- .../edu/sc/seis/TauP/cmdline/TauP_Curve.java | 31 +--- .../edu/sc/seis/TauP/cmdline/TauP_Find.java | 170 +++++++++++------- .../edu/sc/seis/TauP/cmdline/TauP_Path.java | 20 +-- .../seis/TauP/cmdline/TauP_PhaseDescribe.java | 7 +- .../edu/sc/seis/TauP/cmdline/TauP_Pierce.java | 11 +- .../sc/seis/TauP/cmdline/TauP_SetMSeed3.java | 45 ++--- .../edu/sc/seis/TauP/cmdline/TauP_SetSac.java | 8 +- .../edu/sc/seis/TauP/cmdline/TauP_Spikes.java | 13 +- .../edu/sc/seis/TauP/cmdline/TauP_Table.java | 22 +-- .../edu/sc/seis/TauP/cmdline/TauP_Time.java | 145 +++++---------- .../sc/seis/TauP/cmdline/TauP_Wavefront.java | 11 +- .../seis/TauP/cmdline/args/DistanceArgs.java | 48 ++++- .../sc/seis/TauP/cmdline/args/ModelArgs.java | 31 +++- .../seis/TauP/cmdline/args/QmlStaxmlArgs.java | 86 +++++++++ .../edu/sc/seis/TauP/SurfaceWaveTest.java | 20 +-- .../java/edu/sc/seis/TauP/TimeTester.java | 8 +- .../edu/sc/seis/TauP/cmdline/ScatterTest.java | 21 ++- .../sc/seis/TauP/cmdline/TimeOptionsTest.java | 3 +- .../seis/TauP/cmdLineTest/taup_curve_--help | 11 +- .../sc/seis/TauP/cmdLineTest/taup_path_--help | 28 +-- .../seis/TauP/cmdLineTest/taup_phase_--help | 9 +- .../seis/TauP/cmdLineTest/taup_pierce_--help | 26 +-- .../TauP/cmdLineTest/taup_refltrans_--help | 7 +- .../seis/TauP/cmdLineTest/taup_setms3_--help | 24 +-- .../seis/TauP/cmdLineTest/taup_setsac_--help | 6 +- .../seis/TauP/cmdLineTest/taup_spikes_--help | 16 +- .../seis/TauP/cmdLineTest/taup_table_--help | 11 +- .../sc/seis/TauP/cmdLineTest/taup_time_--help | 14 +- .../taup_time_-h_10_-p_P_--deg_35_--json | 3 +- ...-h_10_-p_ttall_--deg_35_--mod_ak135_--json | 35 +++- .../TauP/cmdLineTest/taup_wavefront_--help | 11 +- 62 files changed, 905 insertions(+), 575 deletions(-) create mode 100644 src/main/java/edu/sc/seis/TauP/cmdline/args/QmlStaxmlArgs.java diff --git a/docs/manual/taup_curve.html b/docs/manual/taup_curve.html index 1b732464..843dfc1e 100644 --- a/docs/manual/taup_curve.html +++ b/docs/manual/taup_curve.html @@ -3051,10 +3051,11 @@

Quick search

[--mapwidth=<mapwidth>] [--mapwidthunit=<mapWidthUnit>] [--mw=<mw>] [-o=<outFile>] [--prop=<toolPropsFilename>] [--reddeg=deg/s] [--redkm=km/s] [--rel=phase] [-x=type] - [-y=type] [--xminmax=x x]... [--yminmax=y y]... [[-h=depth] - [--mod=<modelname>] [--stadepth=depth] [--scat=s s]...] - [[-p=phase[,phase...]]... [--phasefile=<phaseFiles>]...] - [--gmt | --json | --svg | --text] + [-y=type] [--xminmax=x x]... [--yminmax=y y]... + [[--mod=<modelname>] [--stadepth=depth] [-h=depth[, + depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]... + [--phasefile=<phaseFiles>]...] [--gmt | --json | --svg | + --text] Plot travel time vs. distance and other curves for seismic phases. Options: @@ -3100,7 +3101,7 @@

Quick search

--ylog Y axis is log --yminmax=y y min and max y axis for plotting Model Args - -h, --evdepth, --sourcedepth=depth + -h, --evdepth, --sourcedepth=depth[,depth...] source depth in km --mod, --model=<modelname> use velocity model "modelname" for calculations. diff --git a/docs/manual/taup_find.html b/docs/manual/taup_find.html index 6d955543..7af96012 100644 --- a/docs/manual/taup_find.html +++ b/docs/manual/taup_find.html @@ -182,8 +182,8 @@

Quick search

--max=<maxActions> [--mw=<mw>] [-o=<outFile>] [--prop=<toolPropsFilename>] [--deg=d[,d...]]... [--rayparamdeg=s/deg [s/deg]]... [--rayparamkm=s/km [s/km]]... - [--time=t [t]]... [--exclude=depth...]... [[-h=depth] - [--mod=<modelname>] [--stadepth=depth] [--scat=s s]...] + [--time=t [t]]... [--exclude=depth...]... [[--mod=<modelname>] + [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]... [--phasefile=<phaseFiles>]...] [--json | --text] Find seismic phases in an earth model. @@ -219,7 +219,7 @@

Quick search

--time=t [t] find arrivals within the given range --verbose enable verbose output Model Args - -h, --evdepth, --sourcedepth=depth + -h, --evdepth, --sourcedepth=depth[,depth...] source depth in km --mod, --model=<modelname> use velocity model "modelname" for calculations. diff --git a/docs/manual/taup_path.html b/docs/manual/taup_path.html index 8f3afb1f..045c7ed2 100644 --- a/docs/manual/taup_path.html +++ b/docs/manual/taup_path.html @@ -2217,17 +2217,18 @@

Quick search

[--color=<coloring>] [--ellipflattening=f] [--mapwidth=<mapwidth>] [--mapwidthunit=<mapWidthUnit>] [--maxpathinc=<maxPathInc>] [-o=<outFile>] - [--prop=<toolPropsFilename>] [--xaxis=type] [--yaxis=type] + [--prop=<toolPropsFilename>] [--qml=<quakemlFilename>] + [--staxml=<stationxmlFilename>] [--xaxis=type] [--yaxis=type] [--degminmax=deg deg]... [--depthminmax=km km]... [--evt=l - l]... [--evt=l l]... [--sta=l l]... [--sta=l l]... [[-h=depth] - [--mod=<modelname>] [--stadepth=depth] [--scat=s s]...] - [[-p=phase[,phase...]]... [--phasefile=<phaseFiles>]...] - [[--deg=d[,d...]]... [--exactdegree=d[,d...]]... - [--exactkilometer=km[,km...]]... [--km=km[,km...]]... - [--rayparamdeg=s/deg[,s/deg...]]... [--rayparamidx=i[, - i...]]... [--rayparamkm=s/km[,s/km...]]... [--rayparamrad=s/rad - [,s/rad...]]... [--takeoff=deg[,deg...]]...]... [--gmt | - --json | --svg | --text] + l]... [--evt=l l]... [--sta=l l]... [--sta=l l]... + [[--mod=<modelname>] [--stadepth=depth] [-h=depth[, + depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]... + [--phasefile=<phaseFiles>]...] [[--deg=d[,d...]]... + [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]... + [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]... + [--rayparamidx=i[,i...]]... [--rayparamkm=s/km[,s/km...]]... + [--rayparamrad=s/rad[,s/rad...]]... [--takeoff=deg[, + deg...]]...]... [--gmt | --json | --svg | --text] Plot ray paths, distance vs. depth, through the earth. Options: @@ -2257,6 +2258,11 @@

Quick search

discontinuities like moho, cmb, iocb --prop=<toolPropsFilename> load defaults from properties file + --qml, --quakeml=<quakemlFilename> + QuakeML file to load for earthquake origins to use + --staxml=<stationxmlFilename> + StationXML file to extract station latitudes and + longitudes from --verbose enable verbose output --withtime include time for each path point, no effect for SVG. @@ -2267,7 +2273,7 @@

Quick search

radius No effect for SVG output. Model Args - -h, --evdepth, --sourcedepth=depth + -h, --evdepth, --sourcedepth=depth[,depth...] source depth in km --mod, --model=<modelname> use velocity model "modelname" for calculations. diff --git a/docs/manual/taup_phase.html b/docs/manual/taup_phase.html index b536c1ed..01b1f877 100644 --- a/docs/manual/taup_phase.html +++ b/docs/manual/taup_phase.html @@ -147,9 +147,10 @@

Quick search

The usage is:

Usage: taup phase [--debug] [--help] [--verbose] [-o=<outFile>]
-                  [--prop=<toolPropsFilename>] [[-h=depth] [--mod=<modelname>]
-                  [--stadepth=depth] [--scat=s s]...] [[-p=phase[,phase...]]...
-                  [--phasefile=<phaseFiles>]...] [--json | --text]
+                  [--prop=<toolPropsFilename>] [[--mod=<modelname>]
+                  [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s s]...]
+                  [[-p=phase[,phase...]]... [--phasefile=<phaseFiles>]...]
+                  [--json | --text]
 Describe a seismic phase in the current model.
 
 Options:
@@ -161,7 +162,7 @@ 
Quick search

                            load defaults from properties file
       --verbose            enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=<modelname>
                            use velocity model "modelname" for calculations.
diff --git a/docs/manual/taup_pierce.html b/docs/manual/taup_pierce.html
index abfff708..2181f4e5 100644
--- a/docs/manual/taup_pierce.html
+++ b/docs/manual/taup_pierce.html
@@ -205,15 +205,16 @@ 
Quick search

                    [--verbose] [--az=<azimuth>] [--az=<azimuth>]
                    [--baz=<backAzimuth>] [--baz=<backAzimuth>]
                    [--ellipflattening=f] [--mw=<mw>] [-o=<outFile>]
-                   [--prop=<toolPropsFilename>] [--pierce=depth[,depth...]]...
-                   [--rel=phase[,phase...]]... [--evt=l l]... [--evt=l l]...
-                   [--sta=l l]... [--sta=l l]... [[-h=depth]
-                   [--mod=<modelname>] [--stadepth=depth] [--scat=s s]...]
-                   [[-p=phase[,phase...]]... [--phasefile=<phaseFiles>]...]
-                   [[--deg=d[,d...]]... [--exactdegree=d[,d...]]...
-                   [--exactkilometer=km[,km...]]... [--km=km[,km...]]...
-                   [--rayparamdeg=s/deg[,s/deg...]]... [--rayparamidx=i[,
-                   i...]]... [--rayparamkm=s/km[,s/km...]]...
+                   [--prop=<toolPropsFilename>] [--qml=<quakemlFilename>]
+                   [--staxml=<stationxmlFilename>] [--pierce=depth[,
+                   depth...]]... [--rel=phase[,phase...]]... [--evt=l l]...
+                   [--evt=l l]... [--sta=l l]... [--sta=l l]...
+                   [[--mod=<modelname>] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                   [--phasefile=<phaseFiles>]...] [[--deg=d[,d...]]...
+                   [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
+                   [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
+                   [--rayparamidx=i[,i...]]... [--rayparamkm=s/km[,s/km...]]...
                    [--rayparamrad=s/rad[,s/rad...]]... [--takeoff=deg[,
                    deg...]]...]... [--json | --text]
 Calculate pierce points for phases at discontinuities in the model.
@@ -241,16 +242,21 @@ 
Quick search

                              additional depth for calculating pierce points
       --prop=<toolPropsFilename>
                              load defaults from properties file
+      --qml, --quakeml=<quakemlFilename>
+                             QuakeML file to load for earthquake origins to use
       --rayp, --onlyrayp     only output the ray parameter
       --rel=phase[,phase...] times relative to the first of the given phases
       --rev                  only prints underside and bottom turn points, e.g.
                                ^ and v
+      --staxml=<stationxmlFilename>
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --time, --onlytime     only output travel time
       --turn                 only prints bottom turning points, e.g. v
       --under                only prints underside reflection points, e.g. ^
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=<modelname>
                              use velocity model "modelname" for calculations.
diff --git a/docs/manual/taup_refltrans.html b/docs/manual/taup_refltrans.html
index d2164f6a..7be8f89d 100644
--- a/docs/manual/taup_refltrans.html
+++ b/docs/manual/taup_refltrans.html
@@ -3311,8 +3311,9 @@ 
Quick search

                       [-o=<outFile>] [--prop=<toolPropsFilename>]
                       [--rpstep=s/km] [-x=type] [--xminmax=x x]... [--yminmax=y
                       y]... [--layer=v v v v v v]... [-y=type...]... [--gmt |
-                      --json | --svg | --text] [[-h=depth] [--mod=<modelname>]
-                      [--stadepth=depth] [--scat=s s]...]
+                      --json | --svg | --text] [[--mod=<modelname>]
+                      [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s
+                      s]...]
 Plot reflection and transmission coefficients for a discontinuity.
 
 Options:
@@ -3360,7 +3361,7 @@ 
Quick search

       --svg                 outputs as SVG
       --text                outputs as Text
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=<modelname>
                             use velocity model "modelname" for calculations.
diff --git a/docs/manual/taup_setms3.html b/docs/manual/taup_setms3.html
index 66887e3f..c26b17ec 100644
--- a/docs/manual/taup_setms3.html
+++ b/docs/manual/taup_setms3.html
@@ -140,13 +140,13 @@ 
Quick search

The usage is:

-
Usage: taup setms3 [--debug] [--geodetic] [--help] [--verbose] [--eh=<ehKey>]
+
Usage: taup setms3 [--debug] [--geodetic] [--help] [--verbose]
                    [--ellipflattening=f] [--prop=<toolPropsFilename>]
                    [--qml=<quakemlFilename>] [--qmltol=<quakeOTimeTol>]
-                   [--staxml=<stationxmlFilename>] [[-h=depth]
-                   [--mod=<modelname>] [--stadepth=depth] [--scat=s s]...]
-                   [[-p=phase[,phase...]]... [--phasefile=<phaseFiles>]...]
-                   [mseed3file...]
+                   [--staxml=<stationxmlFilename>] [--taupeh=<ehKey>]
+                   [[--mod=<modelname>] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                   [--phasefile=<phaseFiles>]...] [mseed3file...]
 Save travel times in the extra header of miniseed3 files.
     https://crotwell.github.io/ms3eh/
 has details on the JSON structure.
@@ -156,9 +156,6 @@ 
Quick search

 Options:
 
       --debug               enable debugging output
-      --eh=<ehKey>          key to store full TauP JSON output within extra
-                              headers within, otherwise use abbreviated 'bag'
-                              markers
       --ellipflattening=f   Elliptical flattening for distance calculations
                               when --geodetic, defaults to WGS84 ~ 1/298.257
       --geodetic            use geodetic latitude for distance calculations,
@@ -170,16 +167,19 @@ 
Quick search

       --prop=<toolPropsFilename>
                             load defaults from properties file
       --qml, --quakeml=<quakemlFilename>
-                            QuakeML file to load to search for origins that
-                              match this waveform
+                            QuakeML file to load for earthquake origins to use
       --qmltol=<quakeOTimeTol>
                             time window to search for origins in a QuakeML file
                               as an ISO8601 string, default value is PT1H.
       --staxml=<stationxmlFilename>
-                            StationXML file to extract station lat/lon from
+                            StationXML file to extract station latitudes and
+                              longitudes from
+      --taupeh=<ehKey>      key to store full TauP JSON output within extra
+                              headers within, otherwise use abbreviated 'bag'
+                              style markers.
       --verbose             enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=<modelname>
                             use velocity model "modelname" for calculations.
diff --git a/docs/manual/taup_setsac.html b/docs/manual/taup_setsac.html
index c5afaba0..d60da8d2 100644
--- a/docs/manual/taup_setsac.html
+++ b/docs/manual/taup_setsac.html
@@ -178,8 +178,8 @@ 
Quick search

 
The usage is:

 
Usage: taup setsac [--debug] [--evdpkm] [--geodetic] [--help] [--verbose]
                    [--ellipflattening=f] [--prop=<toolPropsFilename>]
-                   [[-h=depth] [--mod=<modelname>] [--stadepth=depth] [--scat=s
-                   s]...] [[-p=phase[,phase...]]...
+                   [[--mod=<modelname>] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                    [--phasefile=<phaseFiles>]...] [sacfile...]
 Set headers to travel times of phases using depth and distance from SAC files.
       [sacfile...]          SAC files to process
@@ -200,7 +200,7 @@ 
Quick search

                             load defaults from properties file
       --verbose             enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=<modelname>
                             use velocity model "modelname" for calculations.
diff --git a/docs/manual/taup_spikes.html b/docs/manual/taup_spikes.html
index a9447027..149220c9 100644
--- a/docs/manual/taup_spikes.html
+++ b/docs/manual/taup_spikes.html
@@ -130,10 +130,11 @@ 
Quick search

 
Usage: taup spikes [--debug] [--geodetic] [--help] [--verbose] [--az=<azimuth>]
                    [--az=<azimuth>] [--baz=<backAzimuth>] [--baz=<backAzimuth>]
                    [--ellipflattening=f] [--mw=<mw>] [-o=<outFile>]
-                   [--otime=<origintime>] [--prop=<toolPropsFilename>] [--evt=l
-                   l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
-                   [[-h=depth] [--mod=<modelname>] [--stadepth=depth] [--scat=s
-                   s]...] [[-p=phase[,phase...]]...
+                   [--otime=<origintime>] [--prop=<toolPropsFilename>]
+                   [--qml=<quakemlFilename>] [--staxml=<stationxmlFilename>]
+                   [--evt=l l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
+                   [[--mod=<modelname>] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                    [--phasefile=<phaseFiles>]...] [[--deg=d[,d...]]...
                    [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
                    [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
@@ -159,9 +160,14 @@ 
Quick search

       --otime=<origintime>   event origin time, as ISO8601
       --prop=<toolPropsFilename>
                              load defaults from properties file
+      --qml, --quakeml=<quakemlFilename>
+                             QuakeML file to load for earthquake origins to use
+      --staxml=<stationxmlFilename>
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=<modelname>
                              use velocity model "modelname" for calculations.
diff --git a/docs/manual/taup_table.html b/docs/manual/taup_table.html
index a74e77f1..73d4db65 100644
--- a/docs/manual/taup_table.html
+++ b/docs/manual/taup_table.html
@@ -147,10 +147,11 @@ 
Quick search

 locsat style travel time files.

 
The usage is:

 
Usage: taup table [--debug] [--help] [--verbose] [--header=<headerFile>]
-                  [-o=<outFile>] [--prop=<toolPropsFilename>] [[-h=depth]
-                  [--mod=<modelname>] [--stadepth=depth] [--scat=s s]...]
-                  [[-p=phase[,phase...]]... [--phasefile=<phaseFiles>]...]
-                  [--csv | --json | --locsat | --text]
+                  [-o=<outFile>] [--prop=<toolPropsFilename>]
+                  [[--mod=<modelname>] [--stadepth=depth] [-h=depth[,
+                  depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                  [--phasefile=<phaseFiles>]...] [--csv | --json | --locsat |
+                  --text]
 Create a table of travel times for a range of depths and distances.
 
 Options:
@@ -165,7 +166,7 @@ 
Quick search

                            load defaults from properties file
       --verbose            enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=<modelname>
                            use velocity model "modelname" for calculations.
diff --git a/docs/manual/taup_time.html b/docs/manual/taup_time.html
index 3ee2e24d..c6a1db66 100644
--- a/docs/manual/taup_time.html
+++ b/docs/manual/taup_time.html
@@ -150,10 +150,11 @@ 
Quick search

                  [--time] [--verbose] [--az=<azimuth>] [--az=<azimuth>]
                  [--baz=<backAzimuth>] [--baz=<backAzimuth>]
                  [--ellipflattening=f] [--mw=<mw>] [-o=<outFile>]
-                 [--prop=<toolPropsFilename>] [--rel=phase[,phase...]]...
+                 [--prop=<toolPropsFilename>] [--qml=<quakemlFilename>]
+                 [--staxml=<stationxmlFilename>] [--rel=phase[,phase...]]...
                  [--evt=l l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
-                 [[-h=depth] [--mod=<modelname>] [--stadepth=depth] [--scat=s
-                 s]...] [[-p=phase[,phase...]]...
+                 [[--mod=<modelname>] [--stadepth=depth] [-h=depth[,
+                 depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                  [--phasefile=<phaseFiles>]...] [[--deg=d[,d...]]...
                  [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
                  [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
@@ -181,12 +182,17 @@ 
Quick search

   -o, --output=<outFile>     output to file, default is stdout.
       --prop=<toolPropsFilename>
                              load defaults from properties file
+      --qml, --quakeml=<quakemlFilename>
+                             QuakeML file to load for earthquake origins to use
       --rayp, --onlyrayp     only output the ray parameter
       --rel=phase[,phase...] times relative to the first of the given phases
+      --staxml=<stationxmlFilename>
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --time, --onlytime     only output travel time
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=<modelname>
                              use velocity model "modelname" for calculations.
diff --git a/docs/manual/taup_wavefront.html b/docs/manual/taup_wavefront.html
index f61986a5..23ec3fe7 100644
--- a/docs/manual/taup_wavefront.html
+++ b/docs/manual/taup_wavefront.html
@@ -8864,10 +8864,11 @@ 
Quick search

                       [--mapwidthunit=<mapWidthUnit>] [-o=<outFile>]
                       [--prop=<toolPropsFilename>] [--timestep=<timeStep>]
                       [--xaxis=type] [--yaxis=type] [--degminmax=deg deg]...
-                      [--depthminmax=km km]... [[-h=depth] [--mod=<modelname>]
-                      [--stadepth=depth] [--scat=s s]...] [[-p=phase[,
-                      phase...]]... [--phasefile=<phaseFiles>]...] [--gmt |
-                      --json | --svg | --text]
+                      [--depthminmax=km km]... [[--mod=<modelname>]
+                      [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s
+                      s]...] [[-p=phase[,phase...]]...
+                      [--phasefile=<phaseFiles>]...] [--gmt | --json | --svg |
+                      --text]
 Plot wavefronts of seismic phases at steps in time.
 
 Options:
@@ -8900,7 +8901,7 @@ 
Quick search

                               radius
                             No effect for SVG output.
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=<modelname>
                             use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/VERSION.json b/src/doc/sphinx/source/cmdLineHelp/VERSION.json
index ca8b5a75..c246daa2 100644
--- a/src/doc/sphinx/source/cmdLineHelp/VERSION.json
+++ b/src/doc/sphinx/source/cmdLineHelp/VERSION.json
@@ -2,11 +2,11 @@
   "version": "3.0.0-SNAPSHOT5",
   "name": "TauP",
   "group": "edu.sc.seis",
-  "date": "2024-08-30T19:40:02Z",
+  "date": "2024-09-06T13:16:23Z",
   "git": {
-    "revision": "1480",
-    "shortsha": "66a316be",
-    "sha": "66a316be88a8a854362bb267b59e13b2af8d41c6",
-    "date": "2024-08-30T15:38:34-04:00"
+    "revision": "1485",
+    "shortsha": "1464a61b",
+    "sha": "1464a61b783a48dd00494775a2349e85f40c0724",
+    "date": "2024-09-04T11:12:57-04:00"
   }
 }
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_curve.usage b/src/doc/sphinx/source/cmdLineHelp/taup_curve.usage
index e87e4e42..1a493d3d 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_curve.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_curve.usage
@@ -3,10 +3,11 @@ Usage: taup curve [--debug] [--help] [--legend] [--verbose] [--xabs] [--xlog]
                   [--mapwidth=] [--mapwidthunit=]
                   [--mw=] [-o=] [--prop=]
                   [--reddeg=deg/s] [--redkm=km/s] [--rel=phase] [-x=type]
-                  [-y=type] [--xminmax=x x]... [--yminmax=y y]... [[-h=depth]
-                  [--mod=] [--stadepth=depth] [--scat=s s]...]
-                  [[-p=phase[,phase...]]... [--phasefile=]...]
-                  [--gmt | --json | --svg | --text]
+                  [-y=type] [--xminmax=x x]... [--yminmax=y y]...
+                  [[--mod=] [--stadepth=depth] [-h=depth[,
+                  depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                  [--phasefile=]...] [--gmt | --json | --svg |
+                  --text]
 Plot travel time vs. distance and other curves for seismic phases.
 
 Options:
@@ -52,7 +53,7 @@ Options:
       --ylog               Y axis is log
       --yminmax=y y        min and max y axis for plotting
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=
                            use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_find.usage b/src/doc/sphinx/source/cmdLineHelp/taup_find.usage
index cbed9a2c..58b310ed 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_find.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_find.usage
@@ -3,8 +3,8 @@ Usage: taup find [--amp] [--debug] [--help] [--onlynameddiscon] [--pwaveonly]
                  --max= [--mw=] [-o=]
                  [--prop=] [--deg=d[,d...]]...
                  [--rayparamdeg=s/deg [s/deg]]... [--rayparamkm=s/km [s/km]]...
-                 [--time=t [t]]... [--exclude=depth...]... [[-h=depth]
-                 [--mod=] [--stadepth=depth] [--scat=s s]...]
+                 [--time=t [t]]... [--exclude=depth...]... [[--mod=]
+                 [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s s]...]
                  [[-p=phase[,phase...]]... [--phasefile=]...]
                  [--json | --text]
 Find seismic phases in an earth model.
@@ -40,7 +40,7 @@ Options:
       --time=t [t]         find arrivals within the given range
       --verbose            enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=
                            use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_path.usage b/src/doc/sphinx/source/cmdLineHelp/taup_path.usage
index 5ff8b49e..647ac784 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_path.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_path.usage
@@ -4,17 +4,18 @@ Usage: taup path [--debug] [--geodetic] [--help] [--onlynameddiscon]
                  [--color=] [--ellipflattening=f]
                  [--mapwidth=] [--mapwidthunit=]
                  [--maxpathinc=] [-o=]
-                 [--prop=] [--xaxis=type] [--yaxis=type]
+                 [--prop=] [--qml=]
+                 [--staxml=] [--xaxis=type] [--yaxis=type]
                  [--degminmax=deg deg]... [--depthminmax=km km]... [--evt=l
-                 l]... [--evt=l l]... [--sta=l l]... [--sta=l l]... [[-h=depth]
-                 [--mod=] [--stadepth=depth] [--scat=s s]...]
-                 [[-p=phase[,phase...]]... [--phasefile=]...]
-                 [[--deg=d[,d...]]... [--exactdegree=d[,d...]]...
-                 [--exactkilometer=km[,km...]]... [--km=km[,km...]]...
-                 [--rayparamdeg=s/deg[,s/deg...]]... [--rayparamidx=i[,
-                 i...]]... [--rayparamkm=s/km[,s/km...]]... [--rayparamrad=s/rad
-                 [,s/rad...]]... [--takeoff=deg[,deg...]]...]... [--gmt |
-                 --json | --svg | --text]
+                 l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
+                 [[--mod=] [--stadepth=depth] [-h=depth[,
+                 depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                 [--phasefile=]...] [[--deg=d[,d...]]...
+                 [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
+                 [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
+                 [--rayparamidx=i[,i...]]... [--rayparamkm=s/km[,s/km...]]...
+                 [--rayparamrad=s/rad[,s/rad...]]... [--takeoff=deg[,
+                 deg...]]...]... [--gmt | --json | --svg | --text]
 Plot ray paths, distance vs. depth, through the earth.
 
 Options:
@@ -44,6 +45,11 @@ Options:
                                discontinuities like moho, cmb, iocb
       --prop=
                              load defaults from properties file
+      --qml, --quakeml=
+                             QuakeML file to load for earthquake origins to use
+      --staxml=
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --verbose              enable verbose output
       --withtime             include time for each path point, no effect for
                                SVG.
@@ -54,7 +60,7 @@ Options:
                                radius
                              No effect for SVG output.
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=
                              use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_phase.usage b/src/doc/sphinx/source/cmdLineHelp/taup_phase.usage
index a8e62d43..436092e5 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_phase.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_phase.usage
@@ -1,7 +1,8 @@
 Usage: taup phase [--debug] [--help] [--verbose] [-o=]
-                  [--prop=] [[-h=depth] [--mod=]
-                  [--stadepth=depth] [--scat=s s]...] [[-p=phase[,phase...]]...
-                  [--phasefile=]...] [--json | --text]
+                  [--prop=] [[--mod=]
+                  [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s s]...]
+                  [[-p=phase[,phase...]]... [--phasefile=]...]
+                  [--json | --text]
 Describe a seismic phase in the current model.
 
 Options:
@@ -13,7 +14,7 @@ Options:
                            load defaults from properties file
       --verbose            enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=
                            use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_pierce.usage b/src/doc/sphinx/source/cmdLineHelp/taup_pierce.usage
index 775bbaed..58c4f77f 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_pierce.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_pierce.usage
@@ -3,15 +3,16 @@ Usage: taup pierce [--amp] [--debug] [--first] [--geodetic] [--help]
                    [--verbose] [--az=] [--az=]
                    [--baz=] [--baz=]
                    [--ellipflattening=f] [--mw=] [-o=]
-                   [--prop=] [--pierce=depth[,depth...]]...
-                   [--rel=phase[,phase...]]... [--evt=l l]... [--evt=l l]...
-                   [--sta=l l]... [--sta=l l]... [[-h=depth]
-                   [--mod=] [--stadepth=depth] [--scat=s s]...]
-                   [[-p=phase[,phase...]]... [--phasefile=]...]
-                   [[--deg=d[,d...]]... [--exactdegree=d[,d...]]...
-                   [--exactkilometer=km[,km...]]... [--km=km[,km...]]...
-                   [--rayparamdeg=s/deg[,s/deg...]]... [--rayparamidx=i[,
-                   i...]]... [--rayparamkm=s/km[,s/km...]]...
+                   [--prop=] [--qml=]
+                   [--staxml=] [--pierce=depth[,
+                   depth...]]... [--rel=phase[,phase...]]... [--evt=l l]...
+                   [--evt=l l]... [--sta=l l]... [--sta=l l]...
+                   [[--mod=] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                   [--phasefile=]...] [[--deg=d[,d...]]...
+                   [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
+                   [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
+                   [--rayparamidx=i[,i...]]... [--rayparamkm=s/km[,s/km...]]...
                    [--rayparamrad=s/rad[,s/rad...]]... [--takeoff=deg[,
                    deg...]]...]... [--json | --text]
 Calculate pierce points for phases at discontinuities in the model.
@@ -39,16 +40,21 @@ Options:
                              additional depth for calculating pierce points
       --prop=
                              load defaults from properties file
+      --qml, --quakeml=
+                             QuakeML file to load for earthquake origins to use
       --rayp, --onlyrayp     only output the ray parameter
       --rel=phase[,phase...] times relative to the first of the given phases
       --rev                  only prints underside and bottom turn points, e.g.
                                ^ and v
+      --staxml=
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --time, --onlytime     only output travel time
       --turn                 only prints bottom turning points, e.g. v
       --under                only prints underside reflection points, e.g. ^
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=
                              use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_refltrans.usage b/src/doc/sphinx/source/cmdLineHelp/taup_refltrans.usage
index be000ca5..35b0697a 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_refltrans.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_refltrans.usage
@@ -5,8 +5,9 @@ Usage: taup refltrans [--abs] [--debug] [--down] [--energyflux] [--fsrf]
                       [-o=] [--prop=]
                       [--rpstep=s/km] [-x=type] [--xminmax=x x]... [--yminmax=y
                       y]... [--layer=v v v v v v]... [-y=type...]... [--gmt |
-                      --json | --svg | --text] [[-h=depth] [--mod=]
-                      [--stadepth=depth] [--scat=s s]...]
+                      --json | --svg | --text] [[--mod=]
+                      [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s
+                      s]...]
 Plot reflection and transmission coefficients for a discontinuity.
 
 Options:
@@ -54,7 +55,7 @@ Output Type:
       --svg                 outputs as SVG
       --text                outputs as Text
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=
                             use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_setms3.usage b/src/doc/sphinx/source/cmdLineHelp/taup_setms3.usage
index b108b7f3..3b68f533 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_setms3.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_setms3.usage
@@ -1,10 +1,10 @@
-Usage: taup setms3 [--debug] [--geodetic] [--help] [--verbose] [--eh=]
+Usage: taup setms3 [--debug] [--geodetic] [--help] [--verbose]
                    [--ellipflattening=f] [--prop=]
                    [--qml=] [--qmltol=]
-                   [--staxml=] [[-h=depth]
-                   [--mod=] [--stadepth=depth] [--scat=s s]...]
-                   [[-p=phase[,phase...]]... [--phasefile=]...]
-                   [mseed3file...]
+                   [--staxml=] [--taupeh=]
+                   [[--mod=] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                   [--phasefile=]...] [mseed3file...]
 Save travel times in the extra header of miniseed3 files.
     https://crotwell.github.io/ms3eh/
 has details on the JSON structure.
@@ -14,9 +14,6 @@ has details on the JSON structure.
 Options:
 
       --debug               enable debugging output
-      --eh=          key to store full TauP JSON output within extra
-                              headers within, otherwise use abbreviated 'bag'
-                              markers
       --ellipflattening=f   Elliptical flattening for distance calculations
                               when --geodetic, defaults to WGS84 ~ 1/298.257
       --geodetic            use geodetic latitude for distance calculations,
@@ -28,16 +25,19 @@ Options:
       --prop=
                             load defaults from properties file
       --qml, --quakeml=
-                            QuakeML file to load to search for origins that
-                              match this waveform
+                            QuakeML file to load for earthquake origins to use
       --qmltol=
                             time window to search for origins in a QuakeML file
                               as an ISO8601 string, default value is PT1H.
       --staxml=
-                            StationXML file to extract station lat/lon from
+                            StationXML file to extract station latitudes and
+                              longitudes from
+      --taupeh=      key to store full TauP JSON output within extra
+                              headers within, otherwise use abbreviated 'bag'
+                              style markers.
       --verbose             enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=
                             use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_setsac.usage b/src/doc/sphinx/source/cmdLineHelp/taup_setsac.usage
index f5cdf0e0..d1876ac7 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_setsac.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_setsac.usage
@@ -1,7 +1,7 @@
 Usage: taup setsac [--debug] [--evdpkm] [--geodetic] [--help] [--verbose]
                    [--ellipflattening=f] [--prop=]
-                   [[-h=depth] [--mod=] [--stadepth=depth] [--scat=s
-                   s]...] [[-p=phase[,phase...]]...
+                   [[--mod=] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                    [--phasefile=]...] [sacfile...]
 Set headers to travel times of phases using depth and distance from SAC files.
       [sacfile...]          SAC files to process
@@ -22,7 +22,7 @@ Options:
                             load defaults from properties file
       --verbose             enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=
                             use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_spikes.usage b/src/doc/sphinx/source/cmdLineHelp/taup_spikes.usage
index 159bb098..17234654 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_spikes.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_spikes.usage
@@ -1,10 +1,11 @@
 Usage: taup spikes [--debug] [--geodetic] [--help] [--verbose] [--az=]
                    [--az=] [--baz=] [--baz=]
                    [--ellipflattening=f] [--mw=] [-o=]
-                   [--otime=] [--prop=] [--evt=l
-                   l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
-                   [[-h=depth] [--mod=] [--stadepth=depth] [--scat=s
-                   s]...] [[-p=phase[,phase...]]...
+                   [--otime=] [--prop=]
+                   [--qml=] [--staxml=]
+                   [--evt=l l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
+                   [[--mod=] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                    [--phasefile=]...] [[--deg=d[,d...]]...
                    [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
                    [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
@@ -30,9 +31,14 @@ Options:
       --otime=   event origin time, as ISO8601
       --prop=
                              load defaults from properties file
+      --qml, --quakeml=
+                             QuakeML file to load for earthquake origins to use
+      --staxml=
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=
                              use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_table.usage b/src/doc/sphinx/source/cmdLineHelp/taup_table.usage
index d4abf174..992e0963 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_table.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_table.usage
@@ -1,8 +1,9 @@
 Usage: taup table [--debug] [--help] [--verbose] [--header=]
-                  [-o=] [--prop=] [[-h=depth]
-                  [--mod=] [--stadepth=depth] [--scat=s s]...]
-                  [[-p=phase[,phase...]]... [--phasefile=]...]
-                  [--csv | --json | --locsat | --text]
+                  [-o=] [--prop=]
+                  [[--mod=] [--stadepth=depth] [-h=depth[,
+                  depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                  [--phasefile=]...] [--csv | --json | --locsat |
+                  --text]
 Create a table of travel times for a range of depths and distances.
 
 Options:
@@ -17,7 +18,7 @@ Options:
                            load defaults from properties file
       --verbose            enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=
                            use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_time.usage b/src/doc/sphinx/source/cmdLineHelp/taup_time.usage
index 581f5db7..c89c4aeb 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_time.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_time.usage
@@ -2,10 +2,11 @@ Usage: taup time [--amp] [--debug] [--first] [--geodetic] [--help] [--rayp]
                  [--time] [--verbose] [--az=] [--az=]
                  [--baz=] [--baz=]
                  [--ellipflattening=f] [--mw=] [-o=]
-                 [--prop=] [--rel=phase[,phase...]]...
+                 [--prop=] [--qml=]
+                 [--staxml=] [--rel=phase[,phase...]]...
                  [--evt=l l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
-                 [[-h=depth] [--mod=] [--stadepth=depth] [--scat=s
-                 s]...] [[-p=phase[,phase...]]...
+                 [[--mod=] [--stadepth=depth] [-h=depth[,
+                 depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                  [--phasefile=]...] [[--deg=d[,d...]]...
                  [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
                  [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
@@ -33,12 +34,17 @@ Options:
   -o, --output=     output to file, default is stdout.
       --prop=
                              load defaults from properties file
+      --qml, --quakeml=
+                             QuakeML file to load for earthquake origins to use
       --rayp, --onlyrayp     only output the ray parameter
       --rel=phase[,phase...] times relative to the first of the given phases
+      --staxml=
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --time, --onlytime     only output travel time
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=
                              use velocity model "modelname" for calculations.
diff --git a/src/doc/sphinx/source/cmdLineHelp/taup_wavefront.usage b/src/doc/sphinx/source/cmdLineHelp/taup_wavefront.usage
index 90e6c43d..d978e4fe 100644
--- a/src/doc/sphinx/source/cmdLineHelp/taup_wavefront.usage
+++ b/src/doc/sphinx/source/cmdLineHelp/taup_wavefront.usage
@@ -4,10 +4,11 @@ Usage: taup wavefront [--debug] [--help] [--legend] [--negdist]
                       [--mapwidthunit=] [-o=]
                       [--prop=] [--timestep=]
                       [--xaxis=type] [--yaxis=type] [--degminmax=deg deg]...
-                      [--depthminmax=km km]... [[-h=depth] [--mod=]
-                      [--stadepth=depth] [--scat=s s]...] [[-p=phase[,
-                      phase...]]... [--phasefile=]...] [--gmt |
-                      --json | --svg | --text]
+                      [--depthminmax=km km]... [[--mod=]
+                      [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s
+                      s]...] [[-p=phase[,phase...]]...
+                      [--phasefile=]...] [--gmt | --json | --svg |
+                      --text]
 Plot wavefronts of seismic phases at steps in time.
 
 Options:
@@ -40,7 +41,7 @@ Options:
                               radius
                             No effect for SVG output.
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=
                             use velocity model "modelname" for calculations.
diff --git a/src/main/java/edu/sc/seis/TauP/Arrival.java b/src/main/java/edu/sc/seis/TauP/Arrival.java
index b7419252..e45b08aa 100644
--- a/src/main/java/edu/sc/seis/TauP/Arrival.java
+++ b/src/main/java/edu/sc/seis/TauP/Arrival.java
@@ -719,6 +719,8 @@ public void writeJSON(PrintWriter pw, String indent, boolean withAmp) {
         String NL = "\n";
         pw.write(indent+"{"+NL);
         String innerIndent = indent+"  ";
+
+        pw.write(innerIndent+JSONWriter.valueToString("sourcedepth")+": "+JSONWriter.valueToString((float)getSourceDepth())+","+NL);
         pw.write(innerIndent+JSONWriter.valueToString("distdeg")+": "+JSONWriter.valueToString((float)getModuloDistDeg())+","+NL);
         pw.write(innerIndent+JSONWriter.valueToString("phase")+": "+JSONWriter.valueToString(getName())+","+NL);
         pw.write(innerIndent+JSONWriter.valueToString("time")+": "+JSONWriter.valueToString((float)getTime())+","+NL);
@@ -795,6 +797,7 @@ public void writeJSON(PrintWriter pw, String indent, boolean withAmp) {
     public JSONObject asJSONObject() {
         JSONObject a = new JSONObject();
         a.put("distdeg", (float)getModuloDistDeg());
+        a.put("sourcedepth", (float)getSourceDepth());
         a.put("phase", getName());
         a.put("time", (float)getTime());
         a.put("rayparam", (float)(Math.PI / 180.0 * getRayParam()));
diff --git a/src/main/java/edu/sc/seis/TauP/RayCalculateable.java b/src/main/java/edu/sc/seis/TauP/RayCalculateable.java
index 0cc1aae1..72beb2bd 100644
--- a/src/main/java/edu/sc/seis/TauP/RayCalculateable.java
+++ b/src/main/java/edu/sc/seis/TauP/RayCalculateable.java
@@ -24,6 +24,20 @@ public void withStationBackAzimuth(Location sta, double backazimuth) {
     public abstract boolean isLatLonable();
     public abstract LatLonable getLatLonable();
 
+    public boolean hasSourceDepth() {
+        return evtLatLon != null && evtLatLon.hasDepth();
+    }
+    public Double getSourceDepth() {
+        return evtLatLon != null ? evtLatLon.getDepthKm() : null;
+    }
+
+    public boolean hasReceiverDepth() {
+        return staLatLon != null && staLatLon.hasDepth();
+    }
+    public Double getReceiverDepth() {
+        return staLatLon != null ? staLatLon.getDepthKm() : null;
+    }
+
 
     protected Location staLatLon = null;
     protected Location evtLatLon = null;
diff --git a/src/main/java/edu/sc/seis/TauP/Scatterer.java b/src/main/java/edu/sc/seis/TauP/Scatterer.java
index 537ebc85..0866369c 100644
--- a/src/main/java/edu/sc/seis/TauP/Scatterer.java
+++ b/src/main/java/edu/sc/seis/TauP/Scatterer.java
@@ -10,6 +10,11 @@ public Scatterer(double depth, FixedHemisphereDistanceRay dist) {
         this.dist = dist;
     }
 
+    public Double getDistanceDegree() {
+        // we know created via ofDegrees, so do not need radius to calc
+        return dist.degrees;
+    }
+
     public final double depth;
     public final FixedHemisphereDistanceRay dist;
 
diff --git a/src/main/java/edu/sc/seis/TauP/XYPlotOutput.java b/src/main/java/edu/sc/seis/TauP/XYPlotOutput.java
index 866665f3..f8ac38fe 100644
--- a/src/main/java/edu/sc/seis/TauP/XYPlotOutput.java
+++ b/src/main/java/edu/sc/seis/TauP/XYPlotOutput.java
@@ -16,6 +16,7 @@
 
 import static edu.sc.seis.TauP.AxisType.*;
 import static edu.sc.seis.TauP.cmdline.TauP_AbstractPhaseTool.baseResultAsJSONObject;
+import static edu.sc.seis.TauP.cmdline.args.ModelArgs.depthsToString;
 
 public class XYPlotOutput {
 
@@ -305,7 +306,7 @@ public void printAsSvg(PrintWriter writer, String toolname, List cmdLine
 
     public String getTitle() {
         if (title == null && modelArgs != null) {
-            return modelArgs.getModelName() + " (h=" + modelArgs.getSourceDepth() + " km)";
+            return modelArgs.getModelName() + " (h=" + depthsToString(modelArgs.getSourceDepth()) + " km)";
         } else {
             return title;
         }
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_AbstractPhaseTool.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_AbstractPhaseTool.java
index 9497c260..86f6ee8b 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_AbstractPhaseTool.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_AbstractPhaseTool.java
@@ -201,13 +201,13 @@ public static String[] splitPhaseNameList(String phaseList) {
     }
 
     public static JSONObject baseResultAsJSONObject(String modelName,
-                                                    double depth,
+                                                    List depth,
                                                     double receiverDepth,
                                                     List phaseNameList) {
         JSONObject out = new JSONObject();
 
         out.put("model", modelName);
-        out.put("sourcedepth", (float) depth);
+        out.put("sourcedepth",  depth);
         out.put("receiverdepth", (float) receiverDepth);
         if (phaseNameList != null  ) {
             JSONArray outPhases = new JSONArray();
@@ -277,8 +277,8 @@ public void clearPhaseNames() {
         phaseNames = new ArrayList<>();
     }
 
-    public void setSourceDepth(double depth) {
-        this.modelArgs.setSourceDepth( depth);
+    public void setSingleSourceDepth(double depth) {
+        this.modelArgs.setSourceDepth( Collections.singletonList(depth));
         toolProps.put("taup.source.depth", Double.toString(depth));
         clearPhases();
     }
@@ -314,8 +314,8 @@ public String getTauModelName() {
     /**
      * Gets depth corrected TauModel.
      */
-    public TauModel getTauModelDepthCorrected() throws TauModelException {
-        return modelArgs.depthCorrected();
+    public TauModel getTauModelDepthCorrected(double sourceDepth) throws TauModelException {
+        return modelArgs.depthCorrected(sourceDepth);
     }
 
     public void setModelName(String modelName) {
@@ -340,11 +340,57 @@ public void clearPhases() {
         phases = null;
     }
 
-    public List getSeismicPhases() throws TauModelException {
+    /**
+     * Calculates the seismic phases using a possibly new or changed tau model for the given source depth.
+     */
+    public List calcSeismicPhases(double sourceDepth) throws TauModelException {
+        List newPhases = new ArrayList<>();
+        TauModel tModDepth = modelArgs.depthCorrected(sourceDepth);
+        for (PhaseName phaseName : parsePhaseNameList()) {
+            String tempPhaseName = phaseName.getName();
+            // didn't find it precomputed, so recalculate
+            try {
+                List calcPhaseList = SeismicPhaseFactory.createSeismicPhases(
+                        phaseName.getName(),
+                        tModDepth,
+                        sourceDepth,
+                        modelArgs.getReceiverDepth(),
+                        modelArgs.getScatterer(),
+                        isDEBUG());
+                newPhases.addAll(calcPhaseList);
+                for (SeismicPhase seismicPhase : newPhases) {
+                    if (isVerbose()) {
+                        Alert.info(seismicPhase.toString());
+                    }
+                }
+            } catch (ScatterArrivalFailException e) {
+                Alert.warning(e.getMessage() + ", skipping this phase");
+                if (isVerbose() || isDEBUG()) {
+                    e.printStackTrace();
+                }
+            } catch (TauModelException e) {
+                Alert.warning("Error with phase=" + tempPhaseName,
+                        e.getMessage() + "\nSkipping this phase");
+                if (isVerbose() || isDEBUG()) {
+                    e.printStackTrace();
+                }
+            } finally {
+                if (isVerbose()) {
+                    Alert.info("-----------------");
+                }
+            }
+        }
+        return newPhases;
+    }
+
+    public List getSeismicPhases() throws TauPException {
         if (phases == null) {
-            recalcPhases();
+            phases = new ArrayList<>();
+            for (Double sourceDepth : modelArgs.getSourceDepth()) {
+                phases.addAll(calcSeismicPhases(sourceDepth));
+            }
         }
-        return Collections.unmodifiableList(phases);
+        return phases;
     }
 
     /**
@@ -408,14 +454,9 @@ public List parsePhaseList(String phaseList) {
         return out;
     }
 
-    /**
-     * recalculates the given phases using a possibly new or changed tau model.
-     * This should not need to be called by outside classes as it is called by
-     * depthCorrect, and calculate.
-     */
-    protected synchronized void recalcPhases() throws TauModelException {
+    public synchronized List recalcPhases(double sourceDepth) throws TauModelException {
         List newPhases = new ArrayList<>();
-        TauModel tModDepth = modelArgs.depthCorrected();
+        TauModel tModDepth = modelArgs.depthCorrected(sourceDepth);
         for (PhaseName phaseName : parsePhaseNameList()) {
             String tempPhaseName = phaseName.getName();
             // didn't find it precomputed, so recalculate
@@ -423,7 +464,7 @@ protected synchronized void recalcPhases() throws TauModelException {
                 List calcPhaseList = SeismicPhaseFactory.createSeismicPhases(
                         phaseName.getName(),
                         tModDepth,
-                        modelArgs.getSourceDepth(),
+                        sourceDepth,
                         modelArgs.getReceiverDepth(),
                         modelArgs.getScatterer(),
                         isDEBUG());
@@ -450,7 +491,7 @@ protected synchronized void recalcPhases() throws TauModelException {
                 }
             }
         }
-        phases = newPhases;
+        return newPhases;
     }
 
 
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_AbstractRayTool.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_AbstractRayTool.java
index 4fb58af3..2e2c4e40 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_AbstractRayTool.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_AbstractRayTool.java
@@ -3,13 +3,17 @@
 import edu.sc.seis.TauP.*;
 import edu.sc.seis.TauP.cmdline.args.AbstractOutputTypeArgs;
 import edu.sc.seis.TauP.cmdline.args.DistanceArgs;
+import org.json.JSONException;
 import org.json.JSONObject;
+import org.json.JSONWriter;
 import picocli.CommandLine;
 
 import java.io.IOException;
 import java.io.PrintWriter;
 import java.io.StringWriter;
+import java.util.HashSet;
 import java.util.List;
+import java.util.Set;
 
 public abstract class TauP_AbstractRayTool extends TauP_AbstractPhaseTool {
 
@@ -20,61 +24,105 @@ public TauP_AbstractRayTool(AbstractOutputTypeArgs outputTypeArgs) {
         super(outputTypeArgs);
     }
 
-    public static String resultAsJSON(String modelName,
-                                      double depth,
-                                      double receiverDepth,
-                                      List phases,
-                                      List arrivals) {
-        String Q = ""+'"';
-        String COMMA = ",";
-        String QCOMMA = Q+COMMA;
-        String COLON = ": "; // plus space
-        String S = "  ";
-        String QC = Q+COLON;
-        String QCQ = QC+Q;
-        String SQ = S+Q;
-        // use cast to float to limit digits printed
-        StringWriter sw = new StringWriter();
-        PrintWriter out = new PrintWriter(sw);
-        out.println("{");
-        out.println(SQ+"model"+QCQ+modelName+QCOMMA);
-        out.println(SQ+"sourcedepth"+QC+(float)depth+COMMA);
-        out.println(SQ+"receiverdepth"+QC+(float)receiverDepth+COMMA);
-        out.print(SQ+"phases"+Q+": [");
-        for(int p=0; p depthList,
+                                 double receiverDepth,
+                                 List phases,
+                                 List arrivals) {
+        TauP_AbstractRayTool.writeJSON(pw, indent, modelName, depthList, receiverDepth, phases, arrivals,  false, 4.0f);
+    }
+
+    public static void writeJSON(PrintWriter pw, String indent,
+                                 String modelName,
+                                 List depthList,
+                                 double receiverDepth,
+                                 List phases,
+                                 List arrivals,
+                                 boolean withAmplitude,
+                                 float Mw) {
+        String innerIndent = indent+"  ";
+        String NL = "\n";
+        pw.write("{"+NL);
+        pw.write(innerIndent+ JSONWriter.valueToString("model")+": "+JSONWriter.valueToString(modelName)+","+NL);
+        pw.write(innerIndent+JSONWriter.valueToString("sourcedepthlist")+": [");
+
+        boolean firstDp = true;
+        for (Double depth : depthList) {
+            pw.write((firstDp ? "" : ", ") + JSONWriter.valueToString(depth.floatValue()));
+            firstDp = false;
+        }
+
+        pw.write("]," +NL);
+        pw.write(innerIndent+JSONWriter.valueToString("receiverdepth")+": "+JSONWriter.valueToString((float)receiverDepth)+","+NL);
+        pw.write(innerIndent+JSONWriter.valueToString("phases")+": [ ");
+        boolean first = true;
+        for (SeismicPhase phase : phases) {
+            if (first) {
+                first = false;
+            } else {
+                pw.write(", ");
             }
+            pw.write(JSONWriter.valueToString(phase.getName()));
         }
-        out.println(" ]"+COMMA);
-        out.println(SQ+"arrivals"+Q+": [");
-        for(int j = 0; j < arrivals.size(); j++) {
-            Arrival currArrival = arrivals.get(j);
-            JSONObject asjson = currArrival.asJSONObject();
-            out.print(asjson.toString(2));
-            if (j != arrivals.size()-1) {
-                out.print(COMMA);
+        pw.write(" ],"+NL);
+        if (withAmplitude) {
+            pw.write(innerIndent+JSONWriter.valueToString("Mw")+": "+JSONWriter.valueToString(Mw)+","+NL);
+        }
+        pw.write(innerIndent+JSONWriter.valueToString("arrivals")+": ["+NL);
+        first = true;
+        for (Arrival arrival : arrivals) {
+            if (first) {
+                first = false;
+            } else {
+                pw.write(","+NL);
+            }
+            try {
+                arrival.writeJSON(pw, innerIndent + "  ", withAmplitude);
+            } catch (JSONException e) {
+                System.err.println("Error in json: "+ arrival);
+                throw e;
             }
-            out.println();
         }
-        out.println(S+"]");
-        out.print("}");
-        return sw.toString();
+        pw.write(NL);
+        pw.write(innerIndent+"]"+NL);
+        pw.write("}"+NL);
     }
 
     public DistanceArgs getDistanceArgs() {
         return this.distanceArgs;
     }
 
+
+    @Override
+    public List getSeismicPhases() throws TauPException {
+        List phases = super.getSeismicPhases();
+        // add any depths that are part of eventLatLon locations so we calculate those phases as well
+        Set knownDepths = new HashSet<>();
+        for (DistanceRay dr : getDistanceArgs().getDistances()) {
+            if (dr.hasSourceDepth()) {
+                knownDepths.add(dr.getSourceDepth());
+            }
+        }
+        for (SeismicPhase phase : phases) {
+            knownDepths.remove(phase.getSourceDepth());
+        }
+        for (Double depth : knownDepths) {
+            phases.addAll(calcSeismicPhases(depth));
+        }
+        return phases;
+    }
+
     @Override
     public void validateArguments() throws TauPException {
         if (modelArgs.getTauModel() == null) {
             throw new TauModelException("Model for '"+this.getTauModelName()+"' is null, unable to calculate.");
         }
-        if (modelArgs.getTauModel().getRadiusOfEarth() < modelArgs.getSourceDepth()) {
-            throw new TauModelException("Source depth of "+modelArgs.getSourceDepth()+" in '"+this.getTauModelName()
-                    +"' is greater than radius of earth, "+modelArgs.getTauModel().getRadiusOfEarth()+", unable to calculate.");
+        for (double sourceDepth : modelArgs.getSourceDepth()) {
+            if (modelArgs.getTauModel().getRadiusOfEarth() < sourceDepth) {
+                throw new TauModelException("Source depth of " + sourceDepth + " in '" + this.getTauModelName()
+                        + "' is greater than radius of earth, " + modelArgs.getTauModel().getRadiusOfEarth() + ", unable to calculate.");
+            }
         }
         if (modelArgs.getTauModel().getRadiusOfEarth() < modelArgs.getReceiverDepth()) {
             throw new TauModelException("Receiver depth of "+modelArgs.getReceiverDepth()+" in '"+this.getTauModelName()
@@ -92,4 +140,5 @@ public void validateArguments() throws TauPException {
 
 
     public abstract void printResult(PrintWriter out, List arrivalList) throws IOException, TauPException;
+
 }
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Curve.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Curve.java
index 23e737bc..e9732c64 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Curve.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Curve.java
@@ -36,31 +36,9 @@ public void init() throws TauPException {
 
     @Override
     public void start() throws IOException, TauPException {
-        double tempDepth;
         PrintWriter writer = outputTypeArgs.createWriter(spec.commandLine().getOut());
-        if(modelArgs.getSourceDepth() != -1 * Double.MAX_VALUE) {
-            /* enough info given on cmd line, so just do one calc. */
-
-            List xy = calculate(xAxisType, yAxisType);
-            printResult(writer, xy);
-        } else {
-            StreamTokenizer tokenIn = new StreamTokenizer(new InputStreamReader(System.in));
-            tokenIn.parseNumbers();
-            tokenIn.wordChars(',', ',');
-            tokenIn.wordChars('_', '_');
-            System.out.print("Enter Depth: ");
-            tokenIn.nextToken();
-            tempDepth = tokenIn.nval;
-            if(tempDepth < 0.0 || tempDepth > getRadiusOfEarth()) {
-                System.out.println("Depth must be >= 0.0 and "
-                        + "<= tMod.getRadiusOfEarth().\ndepth = " + tempDepth);
-                return;
-            }
-            setSourceDepth(tempDepth);
-
-            List xy = calculate(xAxisType, yAxisType);
-            printResult(writer, xy);
-        }
+        List xy = calculate(xAxisType, yAxisType);
+        printResult(writer, xy);
         writer.close();
     }
 
@@ -78,12 +56,11 @@ public List calculate(AxisType xAxisType, AxisType yAxisType) th
         return xy;
     }
 
-    public List calculateLinear(AxisType xAxisType, AxisType yAxisType) throws TauModelException, VelocityModelException, SlownessModelException {
-        modelArgs.depthCorrected();
+    public List calculateLinear(AxisType xAxisType, AxisType yAxisType) throws TauPException {
         List phaseList = getSeismicPhases();
         List out = new ArrayList<>();
         for (SeismicPhase phase: phaseList) {
-            String phaseDesc = xAxisType+"/"+yAxisType+" "+phase.getName()+" for a source depth of "+modelArgs.getSourceDepth()+" kilometers in the "+modelArgs.getModelName()+" model";
+            String phaseDesc = xAxisType+"/"+yAxisType+" "+phase.getName()+" for a source depth of "+phase.getSourceDepth()+" kilometers in the "+phase.getTauModel().getModelName()+" model";
             String p_or_s = ColoringArgs.BOTH_PSWAVE;
             if (phase.isAllSWave()) {
                 p_or_s = ColoringArgs.SWAVE;
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Find.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Find.java
index 4f90bae9..739834eb 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Find.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Find.java
@@ -68,81 +68,125 @@ public SeismicPhaseWalk createWalker(TauModel tMod, List excludeDepths)
 
     @Override
     public void start() throws IOException, TauPException {
-        TauModel tMod = modelArgs.depthCorrected();
-        List excludeDepths = getExcludedDepths(tMod);
-        List actualExcludeDepths = matchDepthToDiscon(excludeDepths, tMod.getVelocityModel(), excludeDepthTol);
-
-        SeismicPhaseWalk walker = createWalker(tMod, actualExcludeDepths);
-        List walk = walker.findEndingPaths(maxActions);
         List givenPhaseNames = extractPhaseNames("");
-
         List distanceValues = distanceArgs.getRayCalculatables();
-        if((!distanceValues.isEmpty())) {
-            List phaseList = new ArrayList<>();
-            List phaseNameList = new ArrayList<>();
-            for (ProtoSeismicPhase proto : walk) {
-                phaseNameList.add(proto.getName());
-                phaseList.add(proto.asSeismicPhase());
-            }
-            if ( ! phaseArgs.isEmpty()) {
-                List givenPhases = parsePhaseNameList();
-                for (PhaseName pn : givenPhases) {
-                    phaseNameList.add(pn.getName());
-                    phaseList.addAll(SeismicPhaseFactory.createSeismicPhases(pn.getName(),
-                            tMod, modelArgs.getSourceDepth(), modelArgs.getReceiverDepth(),
-                            modelArgs.getScatterer(), isDEBUG()));
-                }
-            }
-            TauP_Time timeTool = new TauP_Time();
-            timeTool.setPhaseNames(phaseNameList);
-            timeTool.modelArgs = modelArgs;
-            timeTool.outputTypeArgs = outputTypeArgs;
-            timeTool.sourceArgs = sourceArgs;
-            timeTool.withAmplitude= withAmplitude;
-            List arrivalList = timeTool.calcAll(phaseList,distanceValues);
-            if (!times.isEmpty()) {
-                double minTime = times.get(0);
-                double maxTime;
-                if (times.size()>1) {
-                    maxTime = times.get(1);
-                } else {
-                    maxTime = minTime + deltaTime;
-                    minTime = minTime - deltaTime;
-                }
-                List timedArrivalList = new ArrayList<>();
-                for (Arrival a : arrivalList) {
-                    if (minTime <= a.getTime() && a.getTime() <= maxTime) {
-                        timedArrivalList.add(a);
-                    }
-                }
-                arrivalList = timedArrivalList;
+        List arrivalList = new ArrayList<>();
+        List allwalk = new ArrayList<>();
+        for (Double sourceDepth : modelArgs.getSourceDepth()) {
+            TauModel tMod = modelArgs.depthCorrected(sourceDepth);
+            List excludeDepths = getExcludedDepths(tMod);
+            List actualExcludeDepths = matchDepthToDiscon(excludeDepths, tMod.getVelocityModel(), excludeDepthTol);
+
+            SeismicPhaseWalk walker = createWalker(tMod, actualExcludeDepths);
+            List walk = walker.findEndingPaths(maxActions);
+
+            if((!distanceValues.isEmpty())) {
+                arrivalList.addAll(findForDist(walk, tMod, distanceValues));
+            } else {
+                allwalk.addAll(findForAllDepth(walk));
             }
+        }
+        if((!distanceValues.isEmpty())) {
             PrintWriter writer = outputTypeArgs.createWriter(spec.commandLine().getOut());
-            timeTool.printResult(writer, arrivalList);
+            printResult(writer, arrivalList);
             writer.flush();
         } else {
-            System.err.println("No distances given, just find phases...");
-            List givenPhases = new ArrayList<>();
-            if ( ! phaseArgs.isEmpty()) {
-                givenPhases = getSeismicPhases();
+            if (outputTypeArgs.isText()) {
+                printResultText(allwalk);
+            } else if (outputTypeArgs.isJSON()) {
+                printResultJson(allwalk);
             }
-            for (SeismicPhase sp : givenPhases) {
-                if (sp instanceof SimpleSeismicPhase) {
-                    SimpleSeismicPhase ssp = (SimpleSeismicPhase) sp;
-                    walk.add(ssp.getProto());
-                } else {
-                    ScatteredSeismicPhase scat = (ScatteredSeismicPhase) sp;
-                    walk.add(scat.getScatteredPhase().getProto());
+        }
+    }
+
+    public List findForDist(List walk, TauModel tMod, List distanceValues) throws IOException, TauPException {
+        List phaseList = new ArrayList<>();
+        List phaseNameList = new ArrayList<>();
+        for (ProtoSeismicPhase proto : walk) {
+            phaseNameList.add(proto.getName());
+            phaseList.add(proto.asSeismicPhase());
+        }
+        if ( ! phaseArgs.isEmpty()) {
+            List givenPhases = parsePhaseNameList();
+            for (PhaseName pn : givenPhases) {
+                phaseNameList.add(pn.getName());
+                phaseList.addAll(SeismicPhaseFactory.createSeismicPhases(pn.getName(),
+                        tMod, tMod.getSourceDepth(), modelArgs.getReceiverDepth(),
+                        modelArgs.getScatterer(), isDEBUG()));
+            }
+        }
+        TauP_Time timeTool = new TauP_Time();
+        timeTool.setPhaseNames(phaseNameList);
+        timeTool.modelArgs = modelArgs;
+        timeTool.outputTypeArgs = outputTypeArgs;
+        timeTool.sourceArgs = sourceArgs;
+        timeTool.withAmplitude= withAmplitude;
+        List arrivalList = timeTool.calcAll(phaseList,distanceValues);
+        if (!times.isEmpty()) {
+            double minTime = times.get(0);
+            double maxTime;
+            if (times.size()>1) {
+                maxTime = times.get(1);
+            } else {
+                maxTime = minTime + deltaTime;
+                minTime = minTime - deltaTime;
+            }
+            List timedArrivalList = new ArrayList<>();
+            for (Arrival a : arrivalList) {
+                if (minTime <= a.getTime() && a.getTime() <= maxTime) {
+                    timedArrivalList.add(a);
                 }
             }
-            if (outputTypeArgs.isText()) {
-                printResultText(walk);
-            } else if (outputTypeArgs.isJSON()) {
-                printResultJson(walk);
+            arrivalList = timedArrivalList;
+        }
+        return arrivalList;
+    }
+
+
+    public void printResult(PrintWriter out, List arrivalList) throws IOException, TauPException {
+        if (getOutputFormat().equals(OutputTypes.JSON)) {
+            TauP_AbstractRayTool.writeJSON(out, "",
+                    getTauModelName(),
+                    modelArgs.getSourceDepth(),
+                    modelArgs.getReceiverDepth(),
+                    getSeismicPhases(),
+                    arrivalList);
+        } else {
+            boolean onlyFirst = false;
+            boolean onlyPrintTime = false;
+            boolean onlyPrintRayP = false;
+            List relativePhaseName = new ArrayList<>();
+            TauP_Time.printArrivalsAsText(out, arrivalList,
+                    modelArgs.getModelName(),
+                    modelArgs.getSourceDepth(),
+                    getReceiverDepth(),
+                    getScatterer(),
+                    onlyFirst, onlyPrintTime, onlyPrintRayP,
+                    isWithAmplitude(), sourceArgs.getMw(),
+                    relativePhaseName);
+        }
+        out.flush();
+    }
+
+    public List findForAllDepth(List walk) throws IOException, TauPException {
+        System.err.println("No distances given, just find phases...");
+        List givenPhases = new ArrayList<>();
+        if ( ! phaseArgs.isEmpty()) {
+            givenPhases = getSeismicPhases();
+        }
+        for (SeismicPhase sp : givenPhases) {
+            if (sp instanceof SimpleSeismicPhase) {
+                SimpleSeismicPhase ssp = (SimpleSeismicPhase) sp;
+                walk.add(ssp.getProto());
+            } else {
+                ScatteredSeismicPhase scat = (ScatteredSeismicPhase) sp;
+                walk.add(scat.getScatteredPhase().getProto());
             }
         }
+        return walk;
     }
 
+
     public List matchDepthToDiscon(List excludeDepth, VelocityModel vMod, double tol) throws NoSuchLayerException {
         List out = new ArrayList<>();
         double[] disconDepths = vMod.getDisconDepths();
@@ -237,7 +281,7 @@ public void destroy() throws TauPException {
 
     @Override
     public void validateArguments() throws TauPException {
-        TauModel tMod = modelArgs.depthCorrected();
+        TauModel tMod = modelArgs.getTauModel();
         List excludeDepths = getExcludedDepths(tMod);
         // throws if cant find depth near discon
         List actualExcludeDepths = matchDepthToDiscon(excludeDepths, tMod.getVelocityModel(), excludeDepthTol);
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Path.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Path.java
index 596f78d1..8a98fb22 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Path.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Path.java
@@ -147,7 +147,6 @@ public String getOutputFormat() {
 	@Override
 	public List calcAll(List phaseList, List shootables) throws TauPException {
 		List arrivals = new ArrayList<>();
-		modelArgs.depthCorrected();
 		for (SeismicPhase phase : phaseList) {
 			for (RayCalculateable shoot : shootables) {
 				arrivals.addAll(shoot.calculate(phase));
@@ -167,7 +166,12 @@ public void destroy() throws TauPException {
 	@Override
 	public void printResult(PrintWriter out, List arrivalList) throws IOException, TauPException {
 		if (getOutputFormat().equals(OutputTypes.JSON)) {
-			printResultJSON(out, arrivalList);
+			TauP_AbstractRayTool.writeJSON(out, "",
+					getTauModelName(),
+					modelArgs.getSourceDepth(),
+					modelArgs.getReceiverDepth(),
+					getSeismicPhases(),
+					arrivalList);
 		} else if (getOutputFormat().equals(OutputTypes.SVG)) {
 			float pixelWidth = (72.0f * getGraphicOutputTypeArgs().mapwidth);
 			printScriptBeginningSVG(out, arrivalList, pixelWidth, distDepthRange, modelArgs, getCmdLineArgs());
@@ -196,7 +200,7 @@ public void printResult(PrintWriter out, List arrivalList) throws IOExc
 			// text/gmt
 			if (getGraphicOutputTypeArgs().isGMT()) {
 				SvgEarth.printGmtScriptBeginning(out, outputTypeArgs.getOutFileBase(),
-						modelArgs.depthCorrected(), outputTypeArgs.mapwidth,
+						modelArgs.getTauModel(), outputTypeArgs.mapwidth,
 						outputTypeArgs.mapWidthUnit, onlyNamedDiscon,
 						toolNameFromClass(this.getClass()), getCmdLineArgs());
 				if (coloring.getColoring() != ColorType.wavetype) {
@@ -283,20 +287,14 @@ public void labelPathsSVG(PrintWriter out, List arrivalList) {
 		out.println("     ");
 	}
 
-	public void printResultJSON(PrintWriter out, List arrivalList) throws TauPException {
-		String s = resultAsJSON(modelArgs.getModelName(), modelArgs.getSourceDepth(), modelArgs.getReceiverDepth(),
-				parsePhaseNameList(), arrivalList);
-		out.println(s);
-	}
-
 	public void printScriptBeginningSVG(PrintWriter out,
 										List arrivalList,
 										float pixelWidth,
 										DistDepthRange distDepthRange,
 										ModelArgs modelArgs,
-										List cmdLineArgs) throws TauModelException {
+										List cmdLineArgs) throws TauPException {
 
-		TauModel tMod = modelArgs.depthCorrected();
+		TauModel tMod = modelArgs.getTauModel();
 		SvgEarthScaling scaleTrans = calcEarthScaleTrans(arrivalList, distDepthRange);
 		String extraCSS = "";
 
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_PhaseDescribe.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_PhaseDescribe.java
index f7d09c0e..1b7ad548 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_PhaseDescribe.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_PhaseDescribe.java
@@ -52,7 +52,6 @@ public String getOutFileExtension() {
     public void start() throws IOException, TauPException {
         PrintWriter writer = outputTypeArgs.createWriter(spec.commandLine().getOut());
         if (!getSeismicPhases().isEmpty()) {
-            modelArgs.depthCorrected();
             printResult(writer);
         } else {
             writer.println("No phases to describe.");
@@ -65,7 +64,7 @@ public void destroy() throws TauPException {
 
     }
 
-    public void printResult(PrintWriter writer) throws TauModelException {
+    public void printResult(PrintWriter writer) throws TauPException {
         if (getOutputFormat().equals(OutputTypes.TEXT)) {
             printResultText(writer);
         } else if (getOutputFormat().equals(OutputTypes.JSON)) {
@@ -76,7 +75,7 @@ public void printResult(PrintWriter writer) throws TauModelException {
         writer.flush();
     }
 
-    public void printResultText(PrintWriter writer) throws TauModelException {
+    public void printResultText(PrintWriter writer) throws TauPException {
         List phaseList = getSeismicPhases();
         for (SeismicPhase phase: phaseList) {
             writer.println(phase.describe());
@@ -95,7 +94,7 @@ public void printResultText(PrintWriter writer) throws TauModelException {
         }
     }
 
-    public void printResultJSON(PrintWriter writer) throws TauModelException {
+    public void printResultJSON(PrintWriter writer) throws TauPException {
         List phaseList = getSeismicPhases();
         writer.println("[");
         boolean first = true;
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Pierce.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Pierce.java
index a2570ed6..45559bab 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Pierce.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Pierce.java
@@ -118,6 +118,10 @@ public List calcAll(List phaseList, List arrivalList) {
+        printPierceAsText(out, arrivalList);
+    }
+
+    public void printPierceAsText(PrintWriter out, List arrivalList) {
         double prevDepth, nextDepth;
         for (Arrival arrival : arrivalList) {
             out.println("> " + arrival.getCommentLine());
@@ -160,13 +164,6 @@ public void printResultText(PrintWriter out, List arrivalList) {
         }
     }
 
-    @Override
-    public void printResultJSON(PrintWriter out, List arrivalList) throws TauPException {
-        String s = resultAsJSON(modelArgs.getModelName(), modelArgs.getSourceDepth(), modelArgs.getReceiverDepth(),
-                parsePhaseNameList(), arrivalList);
-        out.println(s);
-    }
-
     /**
      * checks to see if the given depth has been "added" as a pierce point.
      */
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_SetMSeed3.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_SetMSeed3.java
index b3eebb18..cd3209b1 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_SetMSeed3.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_SetMSeed3.java
@@ -3,6 +3,7 @@
 import edu.sc.seis.TauP.*;
 import edu.sc.seis.TauP.cmdline.args.GeodeticArgs;
 import edu.sc.seis.TauP.cmdline.args.OutputTypes;
+import edu.sc.seis.TauP.cmdline.args.QmlStaxmlArgs;
 import edu.sc.seis.seisFile.Location;
 import edu.sc.seis.seisFile.SeisFileException;
 import edu.sc.seis.seisFile.fdsnws.quakeml.Event;
@@ -55,22 +56,22 @@ public void start() throws IOException, TauPException {
         }
 
         try {
-            if (stationxmlFilename != null ) {
-                FDSNStationXML staxml = FDSNStationXML.loadStationXML(stationxmlFilename);
+            if (qmlStaxmlArgs.hasStationXML() ) {
+                FDSNStationXML staxml = FDSNStationXML.loadStationXML(qmlStaxmlArgs.getStationxmlFilename());
                 networks = staxml.extractAllNetworks();
             }
         } catch (XMLStreamException | SeisFileException e) {
-            throw new TauPException("Unable to process stationxml from "+stationxmlFilename, e);
+            throw new TauPException("Unable to process stationxml from "+qmlStaxmlArgs.getStationxmlFilename(), e);
         }
         try {
-            if (quakemlFilename != null) {
-                FileReader reader = new FileReader(quakemlFilename);
+            if (qmlStaxmlArgs.hasQml()) {
+                FileReader reader = new FileReader(qmlStaxmlArgs.getQuakemlFilename());
                 Quakeml quakeml = Quakeml.loadQuakeML(reader);
                 quakes = quakeml.extractAllEvents();
                 reader.close();
             }
         } catch (XMLStreamException | SeisFileException e) {
-            throw new TauPException("Unable to process quakeml from "+quakemlFilename, e);
+            throw new TauPException("Unable to process quakeml from "+qmlStaxmlArgs.getQuakemlFilename(), e);
         }
         for (String filename : mseed3FileNames) {
             try {
@@ -166,14 +167,14 @@ public void processRecord(MSeed3Record dr3) throws TauPException {
         } else {
             modelArgs.setReceiverDepth(0);
         }
+        double depth = 0;
         if (evLoc != null) {
-            modelArgs.setSourceDepth(evLoc.getDepthKm());
-        } else {
-            modelArgs.setSourceDepth(0);
+            depth = evLoc.getDepthKm();
         }
+        List seismicPhaseList = calcSeismicPhases( depth);
 
         List arrivals = new ArrayList<>();
-        for (SeismicPhase phase : getSeismicPhases()) {
+        for (SeismicPhase phase : seismicPhaseList) {
             arrivals.addAll(rayCalculateable.calculate(phase));
         }
 
@@ -256,25 +257,6 @@ public void setQuakeOTimeTol(Duration quakeOTimeTol) {
         this.quakeOTimeTol = quakeOTimeTol;
     }
 
-    public String getQuakemlFilename() {
-        return quakemlFilename;
-    }
-
-    @CommandLine.Option(names = {"--qml", "--quakeml"},
-            description = "QuakeML file to load to search for origins that match this waveform")
-    public void setQuakemlFilename(String quakemlFilename) {
-        this.quakemlFilename = quakemlFilename;
-    }
-
-    public String getStationxmlFilename() {
-        return stationxmlFilename;
-    }
-
-    @CommandLine.Option(names = "--staxml", description = "StationXML file to extract station lat/lon from")
-    public void setStationxmlFilename(String stationxmlFilename) {
-        this.stationxmlFilename = stationxmlFilename;
-    }
-
     public List getMseed3FileNames() {
         return mseed3FileNames;
     }
@@ -304,8 +286,9 @@ public void setQuakes(List quakes) {
     protected String ehKey = null;
 
     protected Duration quakeOTimeTol = Duration.ofSeconds(3600);
-    protected String quakemlFilename = null;
-    protected String stationxmlFilename = null;
+
+    @CommandLine.Mixin
+    protected QmlStaxmlArgs qmlStaxmlArgs = new QmlStaxmlArgs();
 
     protected List mseed3FileNames = new ArrayList<>();
 
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_SetSac.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_SetSac.java
index 0737f08b..131161b9 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_SetSac.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_SetSac.java
@@ -186,7 +186,7 @@ public void processSacFile(File f) throws IOException, TauPException {
         }
     }
 
-    public void processSacTimeSeries(SacTimeSeries sacFile, String filenameForError) throws SetSacException, TauModelException, SlownessModelException, NoSuchLayerException {
+    public void processSacTimeSeries(SacTimeSeries sacFile, String filenameForError) throws TauPException {
         SacHeader header = sacFile.getHeader();
         if(SacConstants.isUndef(header.getEvdp())) {
             throw new SetSacException("Depth not set in "
@@ -230,7 +230,7 @@ public void processSacTimeSeries(SacTimeSeries sacFile, String filenameForError)
             /* can't get a distance, skipping */
             throw new SetSacException("Can't get a distance, all distance fields are undef in "+filenameForError);
         }
-        if(!((evdpkm && modelArgs.getSourceDepth() == header.getEvdp()) || (!evdpkm && modelArgs.getSourceDepth() == 1000 * header.getEvdp()))) {
+        if(modelArgs.getSourceDepth().size()!= 1 || !((evdpkm && modelArgs.getSourceDepth().get(0) == header.getEvdp()) || (!evdpkm && modelArgs.getSourceDepth().get(0) == 1000 * header.getEvdp()))) {
             if(!evdpkm && header.getEvdp() != 0 && header.getEvdp() < 1000.0) {
                 Alert.warning("Sac header evdp is < 1000 in "
                                       + filenameForError,
@@ -238,9 +238,9 @@ public void processSacTimeSeries(SacTimeSeries sacFile, String filenameForError)
                                       + "(default), you should use the -evdpkm flag");
             }
             if(evdpkm) {
-                setSourceDepth(header.getEvdp());
+                setSingleSourceDepth(header.getEvdp());
             } else {
-                setSourceDepth(header.getEvdp() / 1000.0);
+                setSingleSourceDepth(header.getEvdp() / 1000.0);
             }
         }
         if(isVerbose()) {
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Spikes.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Spikes.java
index 6acf5401..689f31f2 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Spikes.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Spikes.java
@@ -91,6 +91,9 @@ public void validateArguments() throws TauPException {
         if (!getOutputFormat().equals(MS3)) {
             throw new CommandLine.ParameterException(spec.commandLine(), "Unsupported Output Format: "+getOutputFormat());
         }
+        if (modelArgs.getSourceDepth().size() > 1) {
+            throw new CommandLine.ParameterException(spec.commandLine(), "Multiple source depths unsupported: "+modelArgs.getSourceDepth().size());
+        }
     }
 
     @Override
@@ -119,7 +122,6 @@ public List calcSpikes(List degreesList) throws TauPE
 
 
         validateArguments();
-        modelArgs.depthCorrected();
         List phaseList = getSeismicPhases();
         List spikeRecords = new ArrayList<>();
 
@@ -185,6 +187,9 @@ public List calcSpikes(List degreesList) throws TauPE
     public MSeed3EH createEH(double degrees, ZonedDateTime startDateTime, List allArrivals) {
         MSeed3EH eh = new MSeed3EH();
 
+        // assume single source
+        double sourceDepth = modelArgs.getSourceDepth().isEmpty() ? 0 : modelArgs.getSourceDepth().get(0);
+        RayCalculateable rayCalc = allArrivals.get(0).getRayCalculateable();
         Float deg = (float) degrees;
         Float az = distanceArgs.hasAzimuth() ? distanceArgs.getAzimuth().floatValue() : null;
         Float baz = distanceArgs.hasBackAzimuth() ? distanceArgs.getBackAzimuth().floatValue() : null;
@@ -192,7 +197,7 @@ public MSeed3EH createEH(double degrees, ZonedDateTime startDateTime, List calcWKBJ(List degreesList) throws TauPException {
         validateArguments();
-        modelArgs.depthCorrected();
         List phaseList = getSeismicPhases();
         List spikeRecords = new ArrayList<>();
         float[][] sourceTerm = effectiveSourceTerm( sourceArgs.getMw(), (float) deltaT,  1000);
@@ -504,7 +508,6 @@ public DataOutputStream getOutputStream() throws IOException {
     @Override
     public List calcAll(List phaseList, List shootables) throws TauPException {
         List arrivals = new ArrayList<>();
-        modelArgs.depthCorrected();
         for (SeismicPhase phase : phaseList) {
             List phaseArrivals = new ArrayList<>();
             for (RayCalculateable shoot : shootables) {
@@ -544,7 +547,7 @@ public String getOutFileExtension() {
         return outputTypeArgs.getOutFileExtension();
     }
 
-    public List getDistances() {
+    public List getDistances() throws TauPException {
         return distanceArgs.getDistances();
     }
 
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Table.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Table.java
index 6fd82a18..698c06e0 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Table.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Table.java
@@ -37,6 +37,7 @@
 import java.io.PrintWriter;
 import java.io.StreamTokenizer;
 import java.util.ArrayList;
+import java.util.Collections;
 import java.util.List;
 
 import static edu.sc.seis.TauP.cmdline.TauP_Tool.ABREV_SYNOPSIS;
@@ -456,9 +457,9 @@ protected void jsonTable(PrintWriter out) throws TauPException, IOException {
             for (double distance : distances) {
                 List shootables = List.of(DistanceRay.ofDegrees(distance));
                 List arrivals = calcAll(phaseList, shootables);
-                TauP_Time.writeJSON(out, "",
+                TauP_AbstractRayTool.writeJSON(out, "",
                         getTauModelName(),
-                        modelArgs.getSourceDepth(),
+                        Collections.singletonList(depth),
                         modelArgs.getReceiverDepth(),
                         phaseList,
                         arrivals);
@@ -469,13 +470,12 @@ protected void jsonTable(PrintWriter out) throws TauPException, IOException {
 
     protected List recalcPhases(List phaseNames, double depth) throws TauModelException {
         List newPhases = new ArrayList<>();
-        modelArgs.setSourceDepth(depth);
-        TauModel tModDepth = modelArgs.depthCorrected();
+        TauModel tModDepth = modelArgs.depthCorrected(depth);
         for (PhaseName phaseName : phaseNames) {
             List calcPhaseList = SeismicPhaseFactory.createSeismicPhases(
                     phaseName.getName(),
                     tModDepth,
-                    modelArgs.getSourceDepth(),
+                    depth,
                     modelArgs.getReceiverDepth(),
                     modelArgs.getScatterer(),
                     isDEBUG());
@@ -497,7 +497,7 @@ protected void csvTable(PrintWriter out) throws TauPException {
                     double moduloDist = currArrival.getModuloDistDeg();
                     String line = modelArgs.getModelName() + sep
                             + Outputs.formatDistance(moduloDist).trim() + sep
-                            + Outputs.formatDepth(modelArgs.getSourceDepth()).trim() + sep
+                            + Outputs.formatDepth(depth).trim() + sep
                             + currArrival.getName().trim() + sep
                             + Outputs.formatTime(currArrival.getTime()).trim() + sep
                             + Outputs.formatRayParam(Math.PI / 180.0 * currArrival.getRayParam()).trim() + sep
@@ -515,14 +515,14 @@ protected void csvTable(PrintWriter out) throws TauPException {
 
     protected void genericTable(PrintWriter out) throws TauPException {
         for (double depth : depths) {
-            modelArgs.setSourceDepth(depth); // depth correct happens in modelArgs
+            List seismicPhaseList = calcSeismicPhases( depth);
             for (double distance : distances) {
-                List arrivals = calcAll(getSeismicPhases(), List.of(DistanceRay.ofDegrees(distance)));
+                List arrivals = calcAll(seismicPhaseList, List.of(DistanceRay.ofDegrees(distance)));
                 for (Arrival currArrival : arrivals) {
                     double moduloDist = currArrival.getModuloDistDeg();
                     out.print(modelArgs.getModelName() + " "
                             + Outputs.formatDistance(moduloDist) + " "
-                            + Outputs.formatDepth(modelArgs.getSourceDepth()) + " ");
+                            + Outputs.formatDepth(depth) + " ");
                     out.print(currArrival.getName());
                     out.print("  "
                             + Outputs.formatTime(currArrival.getTime())
@@ -571,10 +571,10 @@ protected void locsatTable(PrintWriter out) throws TauPException {
             out.println();
         }
         for (double depth : depths) {
-            modelArgs.setSourceDepth(depth); // depth correct happens in modelArgs
+            List seismicPhaseList = calcSeismicPhases( depth);
             out.println("#  Travel time for z =    " + depth);
             for (double distance : distances) {
-                List arrivals = calcAll(getSeismicPhases(), List.of(DistanceRay.ofDegrees(distance)));
+                List arrivals = calcAll(seismicPhaseList, List.of(DistanceRay.ofDegrees(distance)));
                 String outString = String.format(float15_4, -1.0) + "    none\n";
                 for (Arrival arrival : arrivals) {
                     if (distance > maxDiff
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Time.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Time.java
index f8e617eb..55e44858 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Time.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Time.java
@@ -21,11 +21,8 @@
 import edu.sc.seis.TauP.cmdline.args.OutputTypes;
 import edu.sc.seis.TauP.cmdline.args.SeismicSourceArgs;
 import edu.sc.seis.TauP.cmdline.args.TextOutputTypeArgs;
-import edu.sc.seis.seisFile.Location;
 import org.json.JSONArray;
-import org.json.JSONException;
 import org.json.JSONObject;
-import org.json.JSONWriter;
 import picocli.CommandLine;
 
 import java.io.*;
@@ -114,11 +111,12 @@ public TauP_Time(String modelName) throws TauModelException {
     @Override
     public List calcAll(List phaseList, List shootables) throws TauPException {
         List arrivals = new ArrayList<>();
-        modelArgs.depthCorrected();
         for (SeismicPhase phase : phaseList) {
             List phaseArrivals = new ArrayList<>();
             for (RayCalculateable shoot : shootables) {
-                phaseArrivals.addAll(shoot.calculate(phase));
+                if (! shoot.hasSourceDepth() || shoot.getSourceDepth() == phase.getSourceDepth()) {
+                    phaseArrivals.addAll(shoot.calculate(phase));
+                }
             }
             if (!onlyFirst) {
                 arrivals.addAll(phaseArrivals);
@@ -130,6 +128,7 @@ public List calcAll(List phaseList, List calcAll(List phaseList, List calcSeismicPhases(double sourceDepth) throws TauModelException {
+        List phaseList = super.calcSeismicPhases(sourceDepth);
         relativePhaseList = new ArrayList<>();
         if (!relativePhaseName.isEmpty()) {
             for (String sName : relativePhaseName) {
                 try {
                     List calcRelPhaseList = SeismicPhaseFactory.createSeismicPhases(
                             sName,
-                            modelArgs.depthCorrected(),
-                            modelArgs.getSourceDepth(),
+                            modelArgs.depthCorrected(sourceDepth),
+                            sourceDepth,
                             modelArgs.getReceiverDepth(),
                             modelArgs.getScatterer(),
                             isDEBUG());
@@ -195,6 +193,7 @@ protected synchronized void recalcPhases() throws TauModelException {
                 }
             }
         }
+        return phaseList;
     }
 
     @Override
@@ -208,6 +207,25 @@ public void printResult(PrintWriter out, List arrivalList) throws IOExc
     }
 
     public void printResultText(PrintWriter out, List arrivalList) {
+        printArrivalsAsText(out, arrivalList,
+                modelArgs.getModelName(),
+                modelArgs.getSourceDepth(),
+                getReceiverDepth(),
+                getScatterer(),
+                onlyFirst, onlyPrintTime, onlyPrintRayP,
+                isWithAmplitude(), sourceArgs.getMw(),
+                relativePhaseName);
+    }
+
+    public static void printArrivalsAsText(PrintWriter out,
+                                           List arrivalList,
+                                           String modelName,
+                                           List sourceDepthList,
+                                           double receiverDepth,
+                                           Scatterer scatterer,
+                                           boolean onlyFirst, boolean onlyPrintTime, boolean onlyPrintRayP,
+                                           boolean withAmplitude, double Mw,
+                                           List relativePhaseName) {
         Arrival currArrival;
         int maxNameLength = 5;
         int maxPuristNameLength = 5;
@@ -226,20 +244,20 @@ public void printResultText(PrintWriter out, List arrivalList) {
         //Format phasePuristFormat = new Format("%-" + maxPuristNameLength + "s");
         String phasePuristFormat = "%-" + maxPuristNameLength + "s";
         if(!(onlyPrintRayP || onlyPrintTime)) {
-            String modelLine =  "\nModel: " + modelArgs.getModelName();
-            if (modelArgs.getReceiverDepth() != 0.0) {
-                modelLine += "  Receiver Depth: "+modelArgs.getReceiverDepth()+" km";
+            String modelLine =  "\nModel: " + modelName;
+            if (receiverDepth != 0.0) {
+                modelLine += "  Receiver Depth: "+receiverDepth+" km";
             }
-            if (getScatterer() != null && getScatterer().dist.getDegrees(getRadiusOfEarth()) != 0.0) {
-                modelLine += "  Scatter Depth: "+ getScattererDepth()+" km Dist: "+ getScatterer().dist.getDegrees(getRadiusOfEarth());
+            if (scatterer != null && scatterer.depth != 0.0) {
+                modelLine += "  Scatter Depth: "+ scatterer.depth+" km Dist: "+ scatterer.getDistanceDegree();
             }
             out.println(modelLine);
             String lineOne = "Distance   Depth   " + String.format(phaseFormat, "Phase")
                     + "   Travel    Ray Param  Takeoff  Incident  Purist   "+String.format(phasePuristFormat, "Purist");
             String lineTwo = "  (deg)     (km)   " + String.format(phaseFormat, "Name ")
                     + "   Time (s)  p (s/deg)   (deg)    (deg)   Distance   "+String.format(phasePuristFormat, "Name");
-            if (isWithAmplitude()) {
-                lineOne += "    Amp  ~"+Outputs.formatDistanceNoPad(sourceArgs.getMw())+" Mw";
+            if (withAmplitude) {
+                lineOne += "    Amp  ~"+Outputs.formatDistanceNoPad(Mw)+" Mw";
                 lineTwo += "  Factor PSv   Sh";
             }
             if (!relativePhaseName.isEmpty()) {
@@ -278,13 +296,13 @@ public void printResultText(PrintWriter out, List arrivalList) {
                     out.print("   * ");
                 }
                 out.print(String.format(phasePuristFormat, currArrival.getPuristName()));
-                if (isWithAmplitude()) {
+                if (withAmplitude) {
                     try {
-                        double ampFactorPSV = currArrival.getAmplitudeFactorPSV(sourceArgs.getMoment());
-                        double ampFactorSH = currArrival.getAmplitudeFactorSH(sourceArgs.getMoment());
+                        double ampFactorPSV = currArrival.getAmplitudeFactorPSV(MomentMagnitude.mw_to_N_m(Mw));
+                        double ampFactorSH = currArrival.getAmplitudeFactorSH(MomentMagnitude.mw_to_N_m(Mw));
                         out.print(" " + Outputs.formatAmpFactor(ampFactorPSV) + " " + Outputs.formatAmpFactor(ampFactorSH));
                     } catch (SlownessModelException | TauModelException | VelocityModelException e) {
-                        throw new RuntimeException("SHould not happen", e);
+                        throw new RuntimeException("Should not happen", e);
                     }
                 }
 
@@ -318,11 +336,7 @@ public void printResultText(PrintWriter out, List arrivalList) {
     }
 
     public void printResultJSON(PrintWriter out, List arrivalList) throws TauModelException, TauPException {
-        writeJSON(out, "", arrivalList);
-    }
-
-    public void writeJSON(PrintWriter pw, String indent, List arrivalList) throws TauModelException {
-        writeJSON(pw, indent,
+        TauP_AbstractRayTool.writeJSON(out, "",
                 getTauModelName(),
                 modelArgs.getSourceDepth(),
                 modelArgs.getReceiverDepth(),
@@ -332,64 +346,8 @@ public void writeJSON(PrintWriter pw, String indent, List arrivalList)
                 sourceArgs.getMw());
     }
 
-    public static void writeJSON(PrintWriter pw, String indent,
-                                 String modelName,
-                                 double depth,
-                                 double receiverDepth,
-                                 List phases,
-                                 List arrivals) {
-        writeJSON(pw, indent, modelName, depth, receiverDepth, phases, arrivals,  false, 4.0f);
-    }
-    public static void writeJSON(PrintWriter pw, String indent,
-                                 String modelName,
-                                 double depth,
-                                 double receiverDepth,
-                                 List phases,
-                                 List arrivals,
-                                 boolean withAmplitude,
-                                 float Mw) {
-        String innerIndent = indent+"  ";
-        String NL = "\n";
-        pw.write("{"+NL);
-        pw.write(innerIndent+JSONWriter.valueToString("model")+": "+JSONWriter.valueToString(modelName)+","+NL);
-        pw.write(innerIndent+JSONWriter.valueToString("sourcedepth")+": "+JSONWriter.valueToString((float)depth)+","+NL);
-        pw.write(innerIndent+JSONWriter.valueToString("receiverdepth")+": "+JSONWriter.valueToString((float)receiverDepth)+","+NL);
-        pw.write(innerIndent+JSONWriter.valueToString("phases")+": [ ");
-        boolean first = true;
-        for (SeismicPhase phase : phases) {
-            if (first) {
-                first = false;
-            } else {
-                pw.write(", ");
-            }
-            pw.write(JSONWriter.valueToString(phase.getName()));
-        }
-        pw.write(" ],"+NL);
-        if (withAmplitude) {
-            pw.write(innerIndent+JSONWriter.valueToString("Mw")+": "+JSONWriter.valueToString(Mw)+","+NL);
-        }
-        pw.write(innerIndent+JSONWriter.valueToString("arrivals")+": ["+NL);
-        first = true;
-        for (Arrival arrival : arrivals) {
-            if (first) {
-                first = false;
-            } else {
-                pw.write(","+NL);
-            }
-            try {
-                arrival.writeJSON(pw, innerIndent + "  ", withAmplitude);
-            } catch (JSONException e) {
-                System.err.println("Erro in json: "+ arrival);
-                throw e;
-            }
-        }
-        pw.write(NL);
-        pw.write(innerIndent+"]"+NL);
-        pw.write("}"+NL);
-    }
-
     public static JSONObject resultAsJSONObject(String modelName,
-                                      double depth,
+                                      List depth,
                                       double receiverDepth,
                                       List phases,
                                       List arrivals) {
@@ -406,27 +364,12 @@ public static JSONObject resultAsJSONObject(String modelName,
      * preforms intialization of the tool. Properties are queried for the
      * default model to load, source depth to use, phases to use, etc.
      */
+    @Override
     public void init() throws TauPException {
         super.init();
     }
 
-    public void printHelp() {
-        Alert.info("Enter:\nh for new depth\nr to recalculate\n"
-                + "p to append phases, \nc to clear phases\n"
-                + "l to list phases\n"
-                + "s for new station lat lon\ne for new event lat lon\n"
-                + "a for new azimuth\nb for new back azimuth\n"
-                + "t for takeoff angle\n"
-                + "m for new model or \nq to quit.\n");
-    }
-
-    public void calcAndPrint(List phaseList, List shootables) throws TauPException, IOException {
-        List arrivalList = calcAll(phaseList, shootables);
-        PrintWriter writer = outputTypeArgs.createWriter(spec.commandLine().getOut());
-        printResult(writer, arrivalList);
-        writer.close();
-    }
-
+    @Override
     public void start() throws IOException, TauPException {
         List distanceValues = getDistanceArgs().getRayCalculatables();
         List arrivalList = calcAll(getSeismicPhases(), distanceValues);
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Wavefront.java b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Wavefront.java
index 9d9391be..aed812c8 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Wavefront.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/TauP_Wavefront.java
@@ -160,10 +160,9 @@ scaleTrans, toolNameFromClass(this.getClass()), getCmdLineArgs(),
             SvgEarth.printSvgEnd(out);
         } else {
             // text/gmt
-
             if (getGraphicOutputTypeArgs().isGMT()) {
                 SvgEarth.printGmtScriptBeginning(out, outputTypeArgs.getOutFileBase(),
-                        modelArgs.depthCorrected(), outputTypeArgs.mapwidth,
+                        modelArgs.getTauModel(), outputTypeArgs.mapwidth,
                         outputTypeArgs.mapWidthUnit, onlyNamedDiscon,
                         toolNameFromClass(this.getClass()), getCmdLineArgs());
                 if (coloring.getColoring() == ColorType.none) {
@@ -211,16 +210,18 @@ scaleTrans, toolNameFromClass(this.getClass()), getCmdLineArgs(),
 
 
 
-    public Map> calcIsochron() throws TauModelException {
+    public Map> calcIsochron() throws TauPException {
         Map> out = new HashMap<>();
 
         double maxTime = 0;
         for (SeismicPhase phase : getSeismicPhases()) {
-            if ( ! phase.hasArrivals()) { continue;}
+            if (!phase.hasArrivals()) {
+                continue;
+            }
             maxTime = Math.max(maxTime, phase.getMaxTime());
             Map> wavefrontPathSegments = calcIsochronSegmentsForPhase(phase, timeStep);
             for (Double timeVal : wavefrontPathSegments.keySet()) {
-                if ( ! out.containsKey(timeVal)) {
+                if (!out.containsKey(timeVal)) {
                     out.put(timeVal, new ArrayList<>());
                 }
                 out.get(timeVal).addAll(wavefrontPathSegments.get(timeVal));
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/args/DistanceArgs.java b/src/main/java/edu/sc/seis/TauP/cmdline/args/DistanceArgs.java
index 14d67df1..975d7442 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/args/DistanceArgs.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/args/DistanceArgs.java
@@ -2,16 +2,21 @@
 
 import edu.sc.seis.TauP.*;
 import edu.sc.seis.seisFile.Location;
+import edu.sc.seis.seisFile.fdsnws.quakeml.Event;
+import edu.sc.seis.seisFile.fdsnws.stationxml.Channel;
+import edu.sc.seis.seisFile.fdsnws.stationxml.Network;
+import edu.sc.seis.seisFile.fdsnws.stationxml.Station;
 import picocli.CommandLine;
 import picocli.CommandLine.ArgGroup;
 import picocli.CommandLine.Option;
 
 import java.util.ArrayList;
 import java.util.List;
+import java.util.Map;
 
 public class DistanceArgs {
 
-    public List getDistances() {
+    public List getDistances() throws TauPException {
         List out = new ArrayList<>();
         for (Double d : distArgs.degreesList) {
             out.add(DistanceRay.ofDegrees(d));
@@ -66,6 +71,37 @@ public List getDistances() {
                 }
             }
         }
+        // load pairs of events and sta/chan
+        Map> networks = qmlStaxmlArgs.loadStationXML();
+        List quakes = qmlStaxmlArgs.loadQuakeML();
+        for (Event evt : quakes) {
+            Location evtLoc = new Location(evt);
+            for (Network net : networks.keySet()) {
+                for (Station sta : networks.get(net)) {
+                    List allChans = new ArrayList<>();
+                    for (Channel chan : sta.getChannelList()) {
+                        Location cLoc = new Location(chan);
+                        boolean found = false;
+                        for (Location prev : allChans) {
+                            if (prev.equals(cLoc)) {
+                                found = true;
+                                break;
+                            }
+                        }
+                        if ( ! found) {
+                            allChans.add(cLoc);
+                        }
+                    }
+                    for (Location chanLoc : allChans) {
+                        if (geodeticArgs.isGeodetic()) {
+                            out.add(DistanceRay.ofGeodeticStationEvent(chanLoc, evtLoc, geodeticArgs.getEllipFlattening()));
+                        } else {
+                            out.add(DistanceRay.ofStationEvent(chanLoc, evtLoc));
+                        }
+                    }
+                }
+            }
+        }
         return out;
     }
 
@@ -223,7 +259,7 @@ public List getTakeoffAngleRays() {
         return rpList;
     }
 
-    public List getRayCalculatables() {
+    public List getRayCalculatables() throws TauPException {
         List out = new ArrayList<>();
         out.addAll(getDistances());
         out.addAll(getRayParamDegRays());
@@ -255,7 +291,10 @@ public boolean hasStationLatLon() {
     }
 
     public void validateArguments() {
-        if (distArgs.allEmpty() && (latLonArgs.eventList.isEmpty() || latLonArgs.stationList.isEmpty())) {
+        if (distArgs.allEmpty()
+                && (latLonArgs.eventList.isEmpty() || latLonArgs.stationList.isEmpty())
+                && ! ( qmlStaxmlArgs.hasQml() && qmlStaxmlArgs.hasStationXML())
+        ) {
             throw new IllegalArgumentException("Must specify at least one distance or station, event.");
         }
         for (Double d : distArgs.takeoffAngle) {
@@ -278,6 +317,9 @@ public void validateArguments() {
     @CommandLine.Mixin
     GeodeticArgs geodeticArgs = new GeodeticArgs();
 
+    @CommandLine.Mixin
+    QmlStaxmlArgs qmlStaxmlArgs = new QmlStaxmlArgs();
+
     public List getStationList() {
         return latLonArgs.stationList;
     }
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/args/ModelArgs.java b/src/main/java/edu/sc/seis/TauP/cmdline/args/ModelArgs.java
index c18d3e31..b3799211 100644
--- a/src/main/java/edu/sc/seis/TauP/cmdline/args/ModelArgs.java
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/args/ModelArgs.java
@@ -19,11 +19,14 @@ public void setModelName(String modelName) {
         modelArgsInner.modelname = modelName;
     }
 
-    public double getSourceDepth() {
+    public List getSourceDepth() {
+        if (modelArgsInner.sourceDepth == null || modelArgsInner.sourceDepth.size() == 0) {
+            return List.of(Double.parseDouble(toolProps.getProperty("taup.source.depth", "0.0")));
+        }
         return modelArgsInner.sourceDepth;
     }
-    public void setSourceDepth(double depth) {
-        modelArgsInner.sourceDepth = depth;
+    public void setSourceDepth(List depths) {
+        modelArgsInner.sourceDepth = depths;
     }
     public double getReceiverDepth() {
         return modelArgsInner.receiverDepth;
@@ -32,6 +35,19 @@ public void setReceiverDepth(double depth) {
         modelArgsInner.receiverDepth = depth;
     }
 
+    public static String depthsToString(List depths) {
+        StringBuffer buf = new StringBuffer();
+        boolean first = true;
+        for (Double d : depths) {
+            if ( ! first) {
+                buf.append(", ");
+            }
+            buf.append(d);
+            first = false;
+        }
+        return buf.toString();
+    }
+
     /**
      * Returns the scatterer if set, null otherwise.
      *
@@ -74,12 +90,12 @@ public List getModelSplitDepths() {
         return modelSplitDepths;
     }
 
-    public TauModel depthCorrected() throws TauModelException {
+    public TauModel depthCorrected(double sourceDepth) throws TauModelException {
         if (tModDepth == null) {
             tModDepth = getTauModel();
         }
-        if (tModDepth.getSourceDepth() != getSourceDepth()) {
-            tModDepth = getTauModel().depthCorrect(getSourceDepth());
+        if (tModDepth.getSourceDepth() != sourceDepth) {
+            tModDepth = getTauModel().depthCorrect(sourceDepth);
             if (!tModDepth.isBranchDepth(getReceiverDepth())) {
                 tModDepth = tModDepth.splitBranch(getReceiverDepth());
             }
@@ -129,8 +145,9 @@ static class ModelArgsInner {
         @CommandLine.Option(names={"-h", "--sourcedepth", "--evdepth"},
                 paramLabel = "depth",
                 defaultValue = "0.0",
+                split=",",
                 description = "source depth in km")
-        double sourceDepth = Double.parseDouble(toolProps.getProperty("taup.source.depth", "0.0"));
+        List sourceDepth = new ArrayList<>();
 
         @CommandLine.Option(names = {"--stadepth", "--receiverdepth"},
                 defaultValue = "0.0",
diff --git a/src/main/java/edu/sc/seis/TauP/cmdline/args/QmlStaxmlArgs.java b/src/main/java/edu/sc/seis/TauP/cmdline/args/QmlStaxmlArgs.java
new file mode 100644
index 00000000..be482498
--- /dev/null
+++ b/src/main/java/edu/sc/seis/TauP/cmdline/args/QmlStaxmlArgs.java
@@ -0,0 +1,86 @@
+package edu.sc.seis.TauP.cmdline.args;
+
+import edu.sc.seis.TauP.DistanceRay;
+import edu.sc.seis.TauP.TauPException;
+import edu.sc.seis.seisFile.Location;
+import edu.sc.seis.seisFile.SeisFileException;
+import edu.sc.seis.seisFile.fdsnws.quakeml.Event;
+import edu.sc.seis.seisFile.fdsnws.quakeml.Quakeml;
+import edu.sc.seis.seisFile.fdsnws.stationxml.Channel;
+import edu.sc.seis.seisFile.fdsnws.stationxml.FDSNStationXML;
+import edu.sc.seis.seisFile.fdsnws.stationxml.Network;
+import edu.sc.seis.seisFile.fdsnws.stationxml.Station;
+import picocli.CommandLine;
+
+import javax.xml.stream.XMLStreamException;
+import java.io.FileReader;
+import java.io.IOException;
+import java.time.Duration;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+public class QmlStaxmlArgs {
+
+    public boolean hasQml() {
+        return getQuakemlFilename() != null;
+    }
+
+    public boolean hasStationXML() {
+        return getStationxmlFilename() != null;
+    }
+
+    public Map> loadStationXML() throws TauPException {
+        Map> networks = new HashMap<>();
+        try {
+            if (stationxmlFilename != null ) {
+                FDSNStationXML staxml = FDSNStationXML.loadStationXML(stationxmlFilename);
+                networks = staxml.extractAllNetworks();
+            }
+        } catch (IOException | XMLStreamException | SeisFileException e) {
+            throw new TauPException("Unable to process stationxml from "+stationxmlFilename, e);
+        }
+        return networks;
+    }
+
+    public List loadQuakeML() throws TauPException {
+        List quakes = new ArrayList<>();
+        try {
+            if (quakemlFilename != null) {
+                FileReader reader = new FileReader(quakemlFilename);
+                Quakeml quakeml = Quakeml.loadQuakeML(reader);
+                quakes = quakeml.extractAllEvents();
+                reader.close();
+            }
+        } catch (IOException | XMLStreamException | SeisFileException e) {
+            throw new TauPException("Unable to process quakeml from "+quakemlFilename, e);
+        }
+        return quakes;
+    }
+
+    public String getQuakemlFilename() {
+        return quakemlFilename;
+    }
+
+    @CommandLine.Option(names = {"--qml", "--quakeml"},
+            required = false,
+            description = "QuakeML file to load for earthquake origins to use")
+    public void setQuakemlFilename(String quakemlFilename) {
+        this.quakemlFilename = quakemlFilename;
+    }
+
+    public String getStationxmlFilename() {
+        return stationxmlFilename;
+    }
+
+    @CommandLine.Option(names = "--staxml",
+            required = false,
+            description = "StationXML file to extract station latitudes and longitudes from")
+    public void setStationxmlFilename(String stationxmlFilename) {
+        this.stationxmlFilename = stationxmlFilename;
+    }
+    
+    protected String quakemlFilename = null;
+    protected String stationxmlFilename = null;
+}
diff --git a/src/test/java/edu/sc/seis/TauP/SurfaceWaveTest.java b/src/test/java/edu/sc/seis/TauP/SurfaceWaveTest.java
index 81a90a63..fbbe124c 100644
--- a/src/test/java/edu/sc/seis/TauP/SurfaceWaveTest.java
+++ b/src/test/java/edu/sc/seis/TauP/SurfaceWaveTest.java
@@ -27,10 +27,10 @@ public void testSurfaceWave() throws TauPException {
         float depth = 23;
         float deg = 30;
         taup.setPhaseNames(List.of(phaseName));
-        taup.setSourceDepth(depth);
+        taup.setSingleSourceDepth(depth);
 
 
-        List phaseList = taup.getSeismicPhases();
+        List phaseList = taup.calcSeismicPhases(depth);
         List distanceValues = List.of(DistanceRay.ofDegrees(deg));
         List arrivals = taup.calcAll(phaseList, distanceValues);
         assertTrue(arrivals.size() > 0,
@@ -58,27 +58,27 @@ public void testDistKm() throws TauPException {
         double phaseVel = 3.2;
         String phaseName = phaseVel+"kmps";
         assertEquals("3.2kmps", phaseName);
-        float depth = 0;
-        double distDeg = distKm/taup.getTauModelDepthCorrected().getRadiusOfEarth()*180.0/Math.PI;
+        float sourceDepth = 0;
+        double distDeg = distKm/taup.getTauModelDepthCorrected(sourceDepth).getRadiusOfEarth()*180.0/Math.PI;
         assertEquals(0.8993, distDeg, 0.0001);
         double longWayDistDeg = 360.0-distDeg;
-        double longWayDistKm = longWayDistDeg*taup.getTauModelDepthCorrected().getRadiusOfEarth()*Math.PI/180.0;
+        double longWayDistKm = longWayDistDeg*taup.getTauModelDepthCorrected(sourceDepth).getRadiusOfEarth()*Math.PI/180.0;
         taup.setPhaseNames(List.of(phaseName));
-        taup.setSourceDepth(depth);
-        List phaseList = taup.getSeismicPhases();
+        taup.setSingleSourceDepth(sourceDepth);
+        List phaseList = taup.calcSeismicPhases(sourceDepth);
         List distanceValues = List.of(DistanceRay.ofDegrees(distDeg));
         List arrivals = taup.calcAll(phaseList, distanceValues);
         assertEquals(2, arrivals.size(),
-                phaseName + " has arrivals for depth " + depth + " at dist " + distDeg);
+                phaseName + " has arrivals for depth " + sourceDepth + " at dist " + distDeg);
         assertTrue(arrivals.get(0).getTime() < arrivals.get(1).getTime());
         assertEquals(distKm/phaseVel,
                 arrivals.get(0).getTime(),
                 0.001f,
-                phaseName + " time for depth " + depth + " at dist " + distDeg+"/"+distKm+" "+arrivals.get(0));
+                phaseName + " time for depth " + sourceDepth + " at dist " + distDeg+"/"+distKm+" "+arrivals.get(0));
         // long way around
         assertEquals(longWayDistKm/phaseVel,
                 arrivals.get(1).getTime(),
                 0.001f,
-                phaseName + " time for depth " + depth + " at dist " + longWayDistDeg+"/"+longWayDistKm+" "+arrivals.get(1));
+                phaseName + " time for depth " + sourceDepth + " at dist " + longWayDistDeg+"/"+longWayDistKm+" "+arrivals.get(1));
     }
 }
diff --git a/src/test/java/edu/sc/seis/TauP/TimeTester.java b/src/test/java/edu/sc/seis/TauP/TimeTester.java
index 6cb182dc..be3cc540 100644
--- a/src/test/java/edu/sc/seis/TauP/TimeTester.java
+++ b/src/test/java/edu/sc/seis/TauP/TimeTester.java
@@ -3,6 +3,7 @@
 import edu.sc.seis.TauP.cmdline.args.ModelArgs;
 
 import java.util.ArrayList;
+import java.util.Collections;
 import java.util.List;
 
 /**
@@ -36,20 +37,21 @@ public void setPhaseNames(List phaseList) {
     List phaseNameList;
 
     public void setSourceDepth(double depth) {
-        modelArgs.setSourceDepth(depth);
+        modelArgs.setSourceDepth(Collections.singletonList(depth));
     }
 
     ModelArgs modelArgs = new ModelArgs();
 
     public List getSeismicPhases() throws TauModelException {
-        TauModel tModDepth = modelArgs.depthCorrected();
+        double sourceDepth = modelArgs.getSourceDepth().get(0);
+        TauModel tModDepth = modelArgs.depthCorrected(sourceDepth);
         boolean isDEBUG = false;
         List newPhases = new ArrayList<>();
         for (String ph : phaseNameList) {
             newPhases.addAll(SeismicPhaseFactory.createSeismicPhases(
                     ph,
                     tModDepth,
-                    modelArgs.getSourceDepth(),
+                    sourceDepth,
                     modelArgs.getReceiverDepth(),
                     modelArgs.getScatterer(),
                     isDEBUG));
diff --git a/src/test/java/edu/sc/seis/TauP/cmdline/ScatterTest.java b/src/test/java/edu/sc/seis/TauP/cmdline/ScatterTest.java
index 06e8910a..961f7443 100644
--- a/src/test/java/edu/sc/seis/TauP/cmdline/ScatterTest.java
+++ b/src/test/java/edu/sc/seis/TauP/cmdline/ScatterTest.java
@@ -148,12 +148,13 @@ public void calcScatterDistDegPKoKPTest() {
     public void isLongWayTest() throws TauPException {
         // PKKP is 234 to 287 deg, so is always long way around
         String modelname = "iasp91";
+        double sourceDepth = 0;
         TauP_Time time = new TauP_Time(modelname);
-        time.setSourceDepth(0);
+        time.setSingleSourceDepth(sourceDepth);
         time.clearPhaseNames();
         time.appendPhaseName("PKKP");
 
-        List arrivals = time.calcAll(time.getSeismicPhases(), List.of(DistanceRay.ofDegrees(70)));
+        List arrivals = time.calcAll(time.calcSeismicPhases(sourceDepth), List.of(DistanceRay.ofDegrees(70)));
         for (Arrival a : arrivals) {
             assertTrue(a.isLongWayAround());
         }
@@ -163,7 +164,8 @@ public void isLongWayTest() throws TauPException {
     public void pierceScatterPKoKP() throws TauPException {
         String modelname = "iasp91";
         TauP_Pierce pierce = new TauP_Pierce(modelname);
-        pierce.setSourceDepth(0);
+        double sourceDepth = 0;
+        pierce.setSingleSourceDepth(sourceDepth);
         Scatterer scat = new Scatterer(3500, 120);
         pierce.setScatterer(scat);
         pierce.clearPhaseNames();
@@ -172,7 +174,7 @@ public void pierceScatterPKoKP() throws TauPException {
 
         DistanceRay dRay = DistanceRay.ofDegrees(360-110);
         ScatteredSeismicPhase scatPhase = (ScatteredSeismicPhase) SeismicPhaseFactory.createSeismicPhases(phaseName,
-                pierce.getTauModelDepthCorrected(),
+                pierce.getTauModelDepthCorrected(sourceDepth),
                 0,
                 0,
                 scat,
@@ -189,17 +191,17 @@ public void pierceScatterPKoKP() throws TauPException {
         List dRayArrivalList = dRay.calcScatteredPhase(scatPhase);
         assertEquals(1, dRayArrivalList.size());
 
-        List arrivalsAt250a = pierce.calcAll(pierce.getSeismicPhases(), List.of(DistanceRay.ofDegrees(360-110)));
+        List arrivalsAt250a = pierce.calcAll(pierce.calcSeismicPhases(sourceDepth), List.of(DistanceRay.ofDegrees(360-110)));
         assertEquals(1, arrivalsAt250a.size());
 
 
-        List arrivals = pierce.calcAll(pierce.getSeismicPhases(), List.of(DistanceRay.ofDegrees(-110)));
+        List arrivals = pierce.calcAll(pierce.calcSeismicPhases(sourceDepth), List.of(DistanceRay.ofDegrees(-110)));
         assertEquals(1, arrivals.size());
         Arrival a_neg110 = arrivals.get(0);
         TimeDist[] p_neg110 = a_neg110.getPierce();
 
 
-        List arrivalsAt250 = pierce.calcAll(pierce.getSeismicPhases(), List.of(DistanceRay.ofDegrees(250)));
+        List arrivalsAt250 = pierce.calcAll(pierce.calcSeismicPhases(sourceDepth), List.of(DistanceRay.ofDegrees(250)));
         assertEquals(1, arrivalsAt250.size());
         Arrival a_250 = arrivalsAt250.get(0);
         TimeDist[] p_250 = a_250.getPierce();
@@ -212,15 +214,16 @@ public void pierceScatterPKoKP() throws TauPException {
     @Test
     public void pathBackscatterPedOP() throws TauPException {
         String modelname = "iasp91";
+        double sourceDepth = 0;
         TauP_Path path = new TauP_Path(modelname);
         Scatterer scat = new Scatterer(800, -10);
         path.setScatterer(scat);
 
         path.clearPhaseNames();
         path.appendPhaseName("PedOP");
-        path.recalcPhases();
+        List seismicPhaseList = path.calcSeismicPhases(sourceDepth);
 
-        List arrivals = path.calcAll(path.getSeismicPhases(), Arrays.asList(DistanceRay.ofDegrees(35)));
+        List arrivals = path.calcAll(seismicPhaseList, Arrays.asList(DistanceRay.ofDegrees(35)));
         assertEquals(1, arrivals.size());
         Arrival a_35 = arrivals.get(0);
         ScatteredArrival scatA = (ScatteredArrival) a_35;
diff --git a/src/test/java/edu/sc/seis/TauP/cmdline/TimeOptionsTest.java b/src/test/java/edu/sc/seis/TauP/cmdline/TimeOptionsTest.java
index 34cbbc72..1bc5d832 100644
--- a/src/test/java/edu/sc/seis/TauP/cmdline/TimeOptionsTest.java
+++ b/src/test/java/edu/sc/seis/TauP/cmdline/TimeOptionsTest.java
@@ -15,11 +15,12 @@ public class TimeOptionsTest {
     @Test
     public void TimeOptions() throws TauPException {
         String modelname = "iasp91";
+        double sourceDepth = 0;
         TauModel tMod = TauModelLoader.load(modelname);
         TauP_Time timeTool = new TauP_Time();
         timeTool.setPhaseNames(List.of( "P" ));
         timeTool.onlyFirst = true;
-        List arrivals = timeTool.calcAll(timeTool.getSeismicPhases(), Collections.singletonList(DistanceRay.ofDegrees(20)));
+        List arrivals = timeTool.calcAll(timeTool.calcSeismicPhases(sourceDepth), Collections.singletonList(DistanceRay.ofDegrees(20)));
 
         assertEquals(1, arrivals.size());
         assertEquals(274.09, arrivals.get(0).getTime(), 0.01);
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_curve_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_curve_--help
index e87e4e42..1a493d3d 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_curve_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_curve_--help
@@ -3,10 +3,11 @@ Usage: taup curve [--debug] [--help] [--legend] [--verbose] [--xabs] [--xlog]
                   [--mapwidth=] [--mapwidthunit=]
                   [--mw=] [-o=] [--prop=]
                   [--reddeg=deg/s] [--redkm=km/s] [--rel=phase] [-x=type]
-                  [-y=type] [--xminmax=x x]... [--yminmax=y y]... [[-h=depth]
-                  [--mod=] [--stadepth=depth] [--scat=s s]...]
-                  [[-p=phase[,phase...]]... [--phasefile=]...]
-                  [--gmt | --json | --svg | --text]
+                  [-y=type] [--xminmax=x x]... [--yminmax=y y]...
+                  [[--mod=] [--stadepth=depth] [-h=depth[,
+                  depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                  [--phasefile=]...] [--gmt | --json | --svg |
+                  --text]
 Plot travel time vs. distance and other curves for seismic phases.
 
 Options:
@@ -52,7 +53,7 @@ Options:
       --ylog               Y axis is log
       --yminmax=y y        min and max y axis for plotting
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=
                            use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_path_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_path_--help
index 5ff8b49e..647ac784 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_path_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_path_--help
@@ -4,17 +4,18 @@ Usage: taup path [--debug] [--geodetic] [--help] [--onlynameddiscon]
                  [--color=] [--ellipflattening=f]
                  [--mapwidth=] [--mapwidthunit=]
                  [--maxpathinc=] [-o=]
-                 [--prop=] [--xaxis=type] [--yaxis=type]
+                 [--prop=] [--qml=]
+                 [--staxml=] [--xaxis=type] [--yaxis=type]
                  [--degminmax=deg deg]... [--depthminmax=km km]... [--evt=l
-                 l]... [--evt=l l]... [--sta=l l]... [--sta=l l]... [[-h=depth]
-                 [--mod=] [--stadepth=depth] [--scat=s s]...]
-                 [[-p=phase[,phase...]]... [--phasefile=]...]
-                 [[--deg=d[,d...]]... [--exactdegree=d[,d...]]...
-                 [--exactkilometer=km[,km...]]... [--km=km[,km...]]...
-                 [--rayparamdeg=s/deg[,s/deg...]]... [--rayparamidx=i[,
-                 i...]]... [--rayparamkm=s/km[,s/km...]]... [--rayparamrad=s/rad
-                 [,s/rad...]]... [--takeoff=deg[,deg...]]...]... [--gmt |
-                 --json | --svg | --text]
+                 l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
+                 [[--mod=] [--stadepth=depth] [-h=depth[,
+                 depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                 [--phasefile=]...] [[--deg=d[,d...]]...
+                 [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
+                 [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
+                 [--rayparamidx=i[,i...]]... [--rayparamkm=s/km[,s/km...]]...
+                 [--rayparamrad=s/rad[,s/rad...]]... [--takeoff=deg[,
+                 deg...]]...]... [--gmt | --json | --svg | --text]
 Plot ray paths, distance vs. depth, through the earth.
 
 Options:
@@ -44,6 +45,11 @@ Options:
                                discontinuities like moho, cmb, iocb
       --prop=
                              load defaults from properties file
+      --qml, --quakeml=
+                             QuakeML file to load for earthquake origins to use
+      --staxml=
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --verbose              enable verbose output
       --withtime             include time for each path point, no effect for
                                SVG.
@@ -54,7 +60,7 @@ Options:
                                radius
                              No effect for SVG output.
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=
                              use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_phase_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_phase_--help
index a8e62d43..436092e5 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_phase_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_phase_--help
@@ -1,7 +1,8 @@
 Usage: taup phase [--debug] [--help] [--verbose] [-o=]
-                  [--prop=] [[-h=depth] [--mod=]
-                  [--stadepth=depth] [--scat=s s]...] [[-p=phase[,phase...]]...
-                  [--phasefile=]...] [--json | --text]
+                  [--prop=] [[--mod=]
+                  [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s s]...]
+                  [[-p=phase[,phase...]]... [--phasefile=]...]
+                  [--json | --text]
 Describe a seismic phase in the current model.
 
 Options:
@@ -13,7 +14,7 @@ Options:
                            load defaults from properties file
       --verbose            enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=
                            use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_pierce_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_pierce_--help
index 775bbaed..58c4f77f 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_pierce_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_pierce_--help
@@ -3,15 +3,16 @@ Usage: taup pierce [--amp] [--debug] [--first] [--geodetic] [--help]
                    [--verbose] [--az=] [--az=]
                    [--baz=] [--baz=]
                    [--ellipflattening=f] [--mw=] [-o=]
-                   [--prop=] [--pierce=depth[,depth...]]...
-                   [--rel=phase[,phase...]]... [--evt=l l]... [--evt=l l]...
-                   [--sta=l l]... [--sta=l l]... [[-h=depth]
-                   [--mod=] [--stadepth=depth] [--scat=s s]...]
-                   [[-p=phase[,phase...]]... [--phasefile=]...]
-                   [[--deg=d[,d...]]... [--exactdegree=d[,d...]]...
-                   [--exactkilometer=km[,km...]]... [--km=km[,km...]]...
-                   [--rayparamdeg=s/deg[,s/deg...]]... [--rayparamidx=i[,
-                   i...]]... [--rayparamkm=s/km[,s/km...]]...
+                   [--prop=] [--qml=]
+                   [--staxml=] [--pierce=depth[,
+                   depth...]]... [--rel=phase[,phase...]]... [--evt=l l]...
+                   [--evt=l l]... [--sta=l l]... [--sta=l l]...
+                   [[--mod=] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                   [--phasefile=]...] [[--deg=d[,d...]]...
+                   [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
+                   [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
+                   [--rayparamidx=i[,i...]]... [--rayparamkm=s/km[,s/km...]]...
                    [--rayparamrad=s/rad[,s/rad...]]... [--takeoff=deg[,
                    deg...]]...]... [--json | --text]
 Calculate pierce points for phases at discontinuities in the model.
@@ -39,16 +40,21 @@ Options:
                              additional depth for calculating pierce points
       --prop=
                              load defaults from properties file
+      --qml, --quakeml=
+                             QuakeML file to load for earthquake origins to use
       --rayp, --onlyrayp     only output the ray parameter
       --rel=phase[,phase...] times relative to the first of the given phases
       --rev                  only prints underside and bottom turn points, e.g.
                                ^ and v
+      --staxml=
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --time, --onlytime     only output travel time
       --turn                 only prints bottom turning points, e.g. v
       --under                only prints underside reflection points, e.g. ^
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=
                              use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_refltrans_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_refltrans_--help
index be000ca5..35b0697a 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_refltrans_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_refltrans_--help
@@ -5,8 +5,9 @@ Usage: taup refltrans [--abs] [--debug] [--down] [--energyflux] [--fsrf]
                       [-o=] [--prop=]
                       [--rpstep=s/km] [-x=type] [--xminmax=x x]... [--yminmax=y
                       y]... [--layer=v v v v v v]... [-y=type...]... [--gmt |
-                      --json | --svg | --text] [[-h=depth] [--mod=]
-                      [--stadepth=depth] [--scat=s s]...]
+                      --json | --svg | --text] [[--mod=]
+                      [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s
+                      s]...]
 Plot reflection and transmission coefficients for a discontinuity.
 
 Options:
@@ -54,7 +55,7 @@ Output Type:
       --svg                 outputs as SVG
       --text                outputs as Text
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=
                             use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_setms3_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_setms3_--help
index b108b7f3..3b68f533 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_setms3_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_setms3_--help
@@ -1,10 +1,10 @@
-Usage: taup setms3 [--debug] [--geodetic] [--help] [--verbose] [--eh=]
+Usage: taup setms3 [--debug] [--geodetic] [--help] [--verbose]
                    [--ellipflattening=f] [--prop=]
                    [--qml=] [--qmltol=]
-                   [--staxml=] [[-h=depth]
-                   [--mod=] [--stadepth=depth] [--scat=s s]...]
-                   [[-p=phase[,phase...]]... [--phasefile=]...]
-                   [mseed3file...]
+                   [--staxml=] [--taupeh=]
+                   [[--mod=] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                   [--phasefile=]...] [mseed3file...]
 Save travel times in the extra header of miniseed3 files.
     https://crotwell.github.io/ms3eh/
 has details on the JSON structure.
@@ -14,9 +14,6 @@ has details on the JSON structure.
 Options:
 
       --debug               enable debugging output
-      --eh=          key to store full TauP JSON output within extra
-                              headers within, otherwise use abbreviated 'bag'
-                              markers
       --ellipflattening=f   Elliptical flattening for distance calculations
                               when --geodetic, defaults to WGS84 ~ 1/298.257
       --geodetic            use geodetic latitude for distance calculations,
@@ -28,16 +25,19 @@ Options:
       --prop=
                             load defaults from properties file
       --qml, --quakeml=
-                            QuakeML file to load to search for origins that
-                              match this waveform
+                            QuakeML file to load for earthquake origins to use
       --qmltol=
                             time window to search for origins in a QuakeML file
                               as an ISO8601 string, default value is PT1H.
       --staxml=
-                            StationXML file to extract station lat/lon from
+                            StationXML file to extract station latitudes and
+                              longitudes from
+      --taupeh=      key to store full TauP JSON output within extra
+                              headers within, otherwise use abbreviated 'bag'
+                              style markers.
       --verbose             enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=
                             use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_setsac_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_setsac_--help
index f5cdf0e0..d1876ac7 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_setsac_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_setsac_--help
@@ -1,7 +1,7 @@
 Usage: taup setsac [--debug] [--evdpkm] [--geodetic] [--help] [--verbose]
                    [--ellipflattening=f] [--prop=]
-                   [[-h=depth] [--mod=] [--stadepth=depth] [--scat=s
-                   s]...] [[-p=phase[,phase...]]...
+                   [[--mod=] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                    [--phasefile=]...] [sacfile...]
 Set headers to travel times of phases using depth and distance from SAC files.
       [sacfile...]          SAC files to process
@@ -22,7 +22,7 @@ Options:
                             load defaults from properties file
       --verbose             enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=
                             use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_spikes_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_spikes_--help
index 159bb098..17234654 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_spikes_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_spikes_--help
@@ -1,10 +1,11 @@
 Usage: taup spikes [--debug] [--geodetic] [--help] [--verbose] [--az=]
                    [--az=] [--baz=] [--baz=]
                    [--ellipflattening=f] [--mw=] [-o=]
-                   [--otime=] [--prop=] [--evt=l
-                   l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
-                   [[-h=depth] [--mod=] [--stadepth=depth] [--scat=s
-                   s]...] [[-p=phase[,phase...]]...
+                   [--otime=] [--prop=]
+                   [--qml=] [--staxml=]
+                   [--evt=l l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
+                   [[--mod=] [--stadepth=depth] [-h=depth[,
+                   depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                    [--phasefile=]...] [[--deg=d[,d...]]...
                    [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
                    [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
@@ -30,9 +31,14 @@ Options:
       --otime=   event origin time, as ISO8601
       --prop=
                              load defaults from properties file
+      --qml, --quakeml=
+                             QuakeML file to load for earthquake origins to use
+      --staxml=
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=
                              use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_table_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_table_--help
index d4abf174..992e0963 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_table_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_table_--help
@@ -1,8 +1,9 @@
 Usage: taup table [--debug] [--help] [--verbose] [--header=]
-                  [-o=] [--prop=] [[-h=depth]
-                  [--mod=] [--stadepth=depth] [--scat=s s]...]
-                  [[-p=phase[,phase...]]... [--phasefile=]...]
-                  [--csv | --json | --locsat | --text]
+                  [-o=] [--prop=]
+                  [[--mod=] [--stadepth=depth] [-h=depth[,
+                  depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
+                  [--phasefile=]...] [--csv | --json | --locsat |
+                  --text]
 Create a table of travel times for a range of depths and distances.
 
 Options:
@@ -17,7 +18,7 @@ Options:
                            load defaults from properties file
       --verbose            enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                            source depth in km
       --mod, --model=
                            use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_--help
index 581f5db7..c89c4aeb 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_--help
@@ -2,10 +2,11 @@ Usage: taup time [--amp] [--debug] [--first] [--geodetic] [--help] [--rayp]
                  [--time] [--verbose] [--az=] [--az=]
                  [--baz=] [--baz=]
                  [--ellipflattening=f] [--mw=] [-o=]
-                 [--prop=] [--rel=phase[,phase...]]...
+                 [--prop=] [--qml=]
+                 [--staxml=] [--rel=phase[,phase...]]...
                  [--evt=l l]... [--evt=l l]... [--sta=l l]... [--sta=l l]...
-                 [[-h=depth] [--mod=] [--stadepth=depth] [--scat=s
-                 s]...] [[-p=phase[,phase...]]...
+                 [[--mod=] [--stadepth=depth] [-h=depth[,
+                 depth...]]... [--scat=s s]...] [[-p=phase[,phase...]]...
                  [--phasefile=]...] [[--deg=d[,d...]]...
                  [--exactdegree=d[,d...]]... [--exactkilometer=km[,km...]]...
                  [--km=km[,km...]]... [--rayparamdeg=s/deg[,s/deg...]]...
@@ -33,12 +34,17 @@ Options:
   -o, --output=     output to file, default is stdout.
       --prop=
                              load defaults from properties file
+      --qml, --quakeml=
+                             QuakeML file to load for earthquake origins to use
       --rayp, --onlyrayp     only output the ray parameter
       --rel=phase[,phase...] times relative to the first of the given phases
+      --staxml=
+                             StationXML file to extract station latitudes and
+                               longitudes from
       --time, --onlytime     only output travel time
       --verbose              enable verbose output
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                              source depth in km
       --mod, --model=
                              use velocity model "modelname" for calculations.
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_-h_10_-p_P_--deg_35_--json b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_-h_10_-p_P_--deg_35_--json
index 938a0d13..34cbd134 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_-h_10_-p_P_--deg_35_--json
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_-h_10_-p_P_--deg_35_--json
@@ -1,10 +1,11 @@
 {
   "model": "iasp91",
-  "sourcedepth": 10,
+  "sourcedepthlist": [10],
   "receiverdepth": 0,
   "phases": [ "P" ],
   "arrivals": [
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "P",
       "time": 412.4329,
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_-h_10_-p_ttall_--deg_35_--mod_ak135_--json b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_-h_10_-p_ttall_--deg_35_--mod_ak135_--json
index 0bc1da5b..cd1c287b 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_-h_10_-p_ttall_--deg_35_--mod_ak135_--json
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_time_-h_10_-p_ttall_--deg_35_--mod_ak135_--json
@@ -1,10 +1,11 @@
 {
   "model": "ak135",
-  "sourcedepth": 10,
+  "sourcedepthlist": [10],
   "receiverdepth": 0,
   "phases": [ "p", "P", "Pn", "Pdiff", "PKP", "PKiKP", "PKIKP", "s", "S", "Sn", "Sdiff", "SKS", "SKIKS", "PcP", "pP", "pPdiff", "pPKP", "pPKIKP", "pPKiKP", "sP", "sPdiff", "sPKP", "sPKIKP", "sPKiKP", "sS", "sSdiff", "sSKS", "sSKIKS", "ScS", "pS", "pSdiff", "pSKS", "pSKIKS", "ScP", "SKP", "SKIKP", "PKKP", "PKIKKIKP", "SKKP", "SKIKKIKP", "PP", "PKPPKP", "PKIKPPKIKP", "SKiKP", "PcS", "PKS", "PKIKS", "PKKS", "PKIKKIKS", "SKKS", "SKIKKIKS", "SKSSKS", "SKIKSSKIKS", "SS", "SP", "PS" ],
   "arrivals": [
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "P",
       "time": 412.51175,
@@ -15,6 +16,7 @@
       "puristname": "P"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "pP",
       "time": 415.59045,
@@ -25,6 +27,7 @@
       "puristname": "pP"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "sP",
       "time": 416.83505,
@@ -35,6 +38,7 @@
       "puristname": "sP"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PP",
       "time": 489.44415,
@@ -45,6 +49,7 @@
       "puristname": "PP"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PP",
       "time": 491.84778,
@@ -55,6 +60,7 @@
       "puristname": "PP"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PP",
       "time": 493.22867,
@@ -65,6 +71,7 @@
       "puristname": "PP"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PP",
       "time": 493.25827,
@@ -75,6 +82,7 @@
       "puristname": "PP"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PP",
       "time": 493.2891,
@@ -85,6 +93,7 @@
       "puristname": "PP"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PcP",
       "time": 564.61035,
@@ -95,6 +104,7 @@
       "puristname": "Pcp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "S",
       "time": 744.42035,
@@ -105,6 +115,7 @@
       "puristname": "S"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "pS",
       "time": 747.98193,
@@ -115,6 +126,7 @@
       "puristname": "pS"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "sS",
       "time": 749.4929,
@@ -125,6 +137,7 @@
       "puristname": "sS"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SP",
       "time": 750.45154,
@@ -135,6 +148,7 @@
       "puristname": "SP"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PS",
       "time": 751.9675,
@@ -145,6 +159,7 @@
       "puristname": "PS"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "ScP",
       "time": 790.18524,
@@ -155,6 +170,7 @@
       "puristname": "Scp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PcS",
       "time": 791.36383,
@@ -165,6 +181,7 @@
       "puristname": "Pcs"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SS",
       "time": 879.71454,
@@ -175,6 +192,7 @@
       "puristname": "SS"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SS",
       "time": 896.2998,
@@ -185,6 +203,7 @@
       "puristname": "SS"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SS",
       "time": 898.70404,
@@ -195,6 +214,7 @@
       "puristname": "SS"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SS",
       "time": 930.7437,
@@ -205,6 +225,7 @@
       "puristname": "SS"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SS",
       "time": 930.7906,
@@ -215,6 +236,7 @@
       "puristname": "SS"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PKiKP",
       "time": 1006.5768,
@@ -225,6 +247,7 @@
       "puristname": "PKikp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "pPKiKP",
       "time": 1010.02234,
@@ -235,6 +258,7 @@
       "puristname": "pPKikp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "sPKiKP",
       "time": 1011.18896,
@@ -245,6 +269,7 @@
       "puristname": "sPKikp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "ScS",
       "time": 1033.835,
@@ -255,6 +280,7 @@
       "puristname": "Scs"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SKiKP",
       "time": 1218.1917,
@@ -265,6 +291,7 @@
       "puristname": "SKikp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PKIKKIKP",
       "time": 1898.0303,
@@ -275,6 +302,7 @@
       "puristname": "PKIkKIkp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SKIKKIKP",
       "time": 2109.5889,
@@ -285,6 +313,7 @@
       "puristname": "SKIkKIkp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PKIKKIKS",
       "time": 2110.7554,
@@ -295,6 +324,7 @@
       "puristname": "PKIkKIks"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SKIKKIKS",
       "time": 2322.2493,
@@ -305,6 +335,7 @@
       "puristname": "SKIkKIks"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PKIKPPKIKP",
       "time": 2405.7,
@@ -315,6 +346,7 @@
       "puristname": "PKIkpPKIkp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "PKPPKP",
       "time": 2501.4019,
@@ -325,6 +357,7 @@
       "puristname": "PKpPKp"
     },
     {
+      "sourcedepth": 10,
       "distdeg": 35,
       "phase": "SKIKSSKIKS",
       "time": 3256.4753,
diff --git a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_wavefront_--help b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_wavefront_--help
index 90e6c43d..d978e4fe 100644
--- a/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_wavefront_--help
+++ b/src/test/resources/edu/sc/seis/TauP/cmdLineTest/taup_wavefront_--help
@@ -4,10 +4,11 @@ Usage: taup wavefront [--debug] [--help] [--legend] [--negdist]
                       [--mapwidthunit=] [-o=]
                       [--prop=] [--timestep=]
                       [--xaxis=type] [--yaxis=type] [--degminmax=deg deg]...
-                      [--depthminmax=km km]... [[-h=depth] [--mod=]
-                      [--stadepth=depth] [--scat=s s]...] [[-p=phase[,
-                      phase...]]... [--phasefile=]...] [--gmt |
-                      --json | --svg | --text]
+                      [--depthminmax=km km]... [[--mod=]
+                      [--stadepth=depth] [-h=depth[,depth...]]... [--scat=s
+                      s]...] [[-p=phase[,phase...]]...
+                      [--phasefile=]...] [--gmt | --json | --svg |
+                      --text]
 Plot wavefronts of seismic phases at steps in time.
 
 Options:
@@ -40,7 +41,7 @@ Options:
                               radius
                             No effect for SVG output.
 Model Args
-  -h, --evdepth, --sourcedepth=depth
+  -h, --evdepth, --sourcedepth=depth[,depth...]
                             source depth in km
       --mod, --model=
                             use velocity model "modelname" for calculations.