diff --git a/Project.toml b/Project.toml
index 69b081d..8d45e28 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,8 +1,9 @@
 name = "JetPackTransforms"
 uuid = "e0630bc5-6e1f-509c-a3e4-0aea24bd2b1c"
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
+DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 JetPack = "24ef3835-3876-54c3-8a7a-956cf69ca0b2"
 Jets = "2a57b368-ab28-5ba9-84aa-637fe7991822"
@@ -10,6 +11,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a"
 
 [compat]
+DSP = "0.8, 0.9"
 FFTW = "1"
 JetPack = "1,2"
 Jets = "1"
diff --git a/src/JetPackTransforms.jl b/src/JetPackTransforms.jl
index 3ce1de4..64dc618 100644
--- a/src/JetPackTransforms.jl
+++ b/src/JetPackTransforms.jl
@@ -1,11 +1,11 @@
 module JetPackTransforms
 
-using FFTW, JetPack, Jets, LinearAlgebra, Wavelets
+using FFTW, JetPack, Jets, LinearAlgebra, Wavelets, DSP, Base.Threads
 
 include("jop_dct.jl")
 include("jop_dwt.jl")
 include("jop_fft.jl")
 include("jop_sft.jl")
-include("jop_slantstack.jl") # requires JopTaper from JetPack
+include("jop_slantstack.jl") # requires JopTaper from JetPack and WaveFD.shiftforward! and WaveFD.shiftadjoint! from WaveFD.jl
 
 end
diff --git a/src/jop_slantstack.jl b/src/jop_slantstack.jl
index 55812db..d1ba160 100644
--- a/src/jop_slantstack.jl
+++ b/src/jop_slantstack.jl
@@ -6,24 +6,23 @@ The domain of the operator is `nz` x `nh` with precision T, `dz` is the depth sp
 `dh` is the offset spacing, and `h0` is the origin of the offset axis.  The additional named optional arguments
 along with their default values are,
 
-* `mode="depth` - choose between "depth" and "time" to specify if the input domain is `z-h` or `t-h`.
-* `theta=-60:1.0:60` - range of incidence angles used when `mode="depth"`.
+* `mode="depth` - choose between "depth" and "time" to specify if the input domain is `z-h` or `t-h`
+* `theta=-45:1.0:45` - range of incidence angles used when `mode="depth"`
 * `p=range(-dz/dh,dz/dh,128)` - ray parameter sampling used when `mode="time"`
 * `padz=0.0,padh=0.0` - fractional padding in depth and offset to apply before applying the Fourier transform
-* `taperz=(0,0)` - beginning and end taper (fractional) in the z-direction (or t-direction) before transforming from `z-h` to `kz-kh`
-* `taperh=(0,0)` - beginning and end taper (fractional) in the h-direction before transforming from `z-h` to `kz-kh`
-* `taperkz=(0,0)` - beginning and end taper (fractional) in the kz-direction (or frequency) before transforming from `kz-kh` to `z-h` or `f-kh` to `t-h`
-* `taperkh=(0,0)` - beginning and end taper (fractional) in the kh-direction before transforming from `kz-kh` to `z-h` or `f-kh` to `t-h`
+* `taperz=(0,0)` - beginning and end taper (fractional) in the z-direction (or t-direction) before transforming from `z-h` to `kz-kh` or `t-h` to `f-kh`
+* `taperh=(0,0)` - beginning and end taper (fractional) in the h-direction before transforming from `z-h` to `kz-kh` or `t-h` to `f-kh`
+* `taperkz=(0,0)` - beginning and end taper (fractional) in the kz-direction (or frequency) before sampling
+* `taperkh=(0,0)` - beginning and end taper (fractional) in the kh-direction before sampling
 
 # Notes
 
-* It should be possible to extend this operator to 3D
 * If the mode is "time", then `padz` is the padding for the time dimension, `dz` is the time sampling interval, `taperz` is the taper for the time dimension and `tapkerkz` is the taper for frequency.
 * For mode="depth", typically theta needs to cover both positive and negative angles and must be < 90 degrees.
 """
 function JopSlantStack(
         dom::JetAbstractSpace{T};
-        theta = collect(-60.0:1.0:60.0),
+        theta = collect(-45.0:1.0:45.0),
         p = nothing,
         dz = 1.0,
         dh = 1.0,
@@ -44,17 +43,11 @@ function JopSlantStack(
 
     # kz
     nzfft = nextprod([2,3,5,7], round(Int, nz*(1 + padz)))
-    kn = pi/dz
-    dk = kn/nzfft
-    kz = dk*[0:div(nzfft,2)+1;]
-
-    # symmetric zero padding for z
-    nz_pad_half = div(nzfft - nz, 2)
-    nz_pad_rng = (nz_pad_half+1):(nz_pad_half+nz)
+    kz = rfftfreq(nzfft, 2*pi / dz)
 
     # kh
     nhfft = nextprod([2,3,5,7], round(Int, nh*(1+padh)))
-    kh = 2 * fftfreq(nhfft, pi/dh)
+    kh = fftfreq(nhfft, 2*pi / dh)
 
     # tan(theta) - used for mode=="depth"
     tant = @. tan(deg2rad(theta))
@@ -71,20 +64,21 @@ function JopSlantStack(
 
     # tapers
     TX = JopTaper(dom, (1,2), (taperz[1],taperh[1]), (taperz[2],taperh[2]))
-    TK = JopTaper(JetSpace(Complex{eltype(dom)},div(nzfft,2)+1,mode=="time" ? length(p) : length(tant)), (1,2), (taperkz[1], taperkh[1]), (taperkz[2], taperkh[2]), mode=(:normal,:fftshift))
+    TK = JopTaper(JetSpace(Complex{eltype(dom)},div(nzfft,2)+1, nhfft), (1,2), (taperkz[1], taperkh[1]), (taperkz[2], taperkh[2]), mode=(:normal,:fftshift))
 
     JopLn(dom = dom, rng = JetSpace(T, nz, mode == "time" ? length(p) : length(tant)), df! = JopSlantStack_df!, df′! = JopSlantStack_df′!,
-        s = (;mode, nzfft, nz_pad_rng, nhfft, kz, kh, tant, p, h0, TX, TK))
+        s = (;mode, nzfft, nhfft, kz, kh, tant, p, h0, TX, TK))
 end
 export JopSlantStack
 
-function JopSlantStack_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; mode, nzfft, nz_pad_rng, nhfft, kz, kh, tant, p, h0, TX, TK, kwargs...) where {T}
+function JopSlantStack_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; mode, nzfft, nhfft, kz, kh, tant, p, h0, TX, TK, kwargs...) where {T}
     nh, np, dh = size(m,2), mode == "time" ? length(p) : length(tant), abs(kh[2]-kh[1])
+    nz = size(d, 1)
 
     mpad = zeros(T, nzfft, nhfft)
-    mpad[nz_pad_rng,1:nh] = TX*m
+    mpad[1:nz,1:nh] = TX*m
 
-    M = rfft(mpad)
+    M = TK*rfft(mpad)
 
     compute_kh = mode == "depth" ? slantstack_compute_kh_from_kz : slantstack_compute_kh_from_frequency
     is_out_of_bounds = (ikh_m1, ikh_p1)->(ikh_m1 < 1 || ikh_p1 > nhfft)
@@ -109,17 +103,18 @@ function JopSlantStack_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; mode, n
         D[ikz,ip] = a_m1*d_m1 + a_p1*d_p1
     end
 
-    _d = brfft(TK*D, nzfft, 1)
-    d .= _d[nz_pad_rng,1:np] ./ nzfft
+    _d = brfft(D, nzfft, 1)
+    d .= _d[1:nz,1:np] ./ nzfft
 end
 
-function JopSlantStack_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; mode, nzfft, nz_pad_rng, nhfft, kz, kh, tant, p, h0, TX, TK, kwargs...) where {T}
+function JopSlantStack_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; mode, nzfft, nhfft, kz, kh, tant, p, h0, TX, TK, kwargs...) where {T}
     nh, np, dh = size(m,2), mode == "time" ? length(p) : length(tant), abs(kh[2]-kh[1])
+    nz = size(d, 1)
 
     dpad = zeros(T, nzfft, np)
-    dpad[nz_pad_rng,:] = d
+    dpad[1:nz,:] = d
 
-    D = TK * (rfft(dpad, 1) ./ nzfft)
+    D = (rfft(dpad, 1) ./ nzfft)
 
     compute_kh = mode == "depth" ? slantstack_compute_kh_from_kz : slantstack_compute_kh_from_frequency
     is_out_of_bounds = (ikh_m1, ikh_p1)->(ikh_m1 < 1 || ikh_p1 > nhfft)
@@ -145,7 +140,7 @@ function JopSlantStack_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; mode
         M[ikz,ikh_m1] += a_m1*m_m1
     end
 
-    m .= TX * (brfft(M, nzfft)[nz_pad_rng,1:nh])
+    m .= TX * (brfft(TK * M, nzfft)[1:nz,1:nh])
 end
 
 @inline function slantstack_compute_kh_from_kz(ikz::Int64, ip::Int64, tant, p, kz, kh, nhfft)
@@ -171,3 +166,791 @@ end
 
     ikh_m1, ikh_p1, _kh
 end
+
+"""
+    A = JopSlantStack3D(dom[; dz=1.0, dhx=1.0, dhy=1.0, hx0=0.0, hy0=0.0, ...])
+
+where `A` is the 3D slant-stack operator mapping for `z-hx-hy` to `z-θ-ϕ` (depth mode) or `t-hx-hy` to `tau-px-py` (time mode).
+The domain of the operator is `nz` x `nhx` x `nhy` with precision T, `dz` is the depth spacing (or time interval),
+`dhx` and `dhy` are the offset spacings, and `hx0` and `hy0` are the origins of the offset axes. The additional named optional arguments
+along with their default values are,
+
+* `mode="depth` - choose between "depth" and "time" to specify if the input domain is `z-hx-hy` or `t-hx-hy`.
+* `theta=-45:1.0:45` - range of incidence angles used when `mode="depth"`.
+* `phi=0:45.0:135` - range of azimuth angles used when `mode="depth"`.
+* `dip=0.0` - geologic dip (degrees) used when `mode="depth"`
+* `azimuth=0.0` - geologic azimuth (degrees) used when `mode="depth"`
+* `px=range(-dz/dhx,dz/dhx,64)` - ray parameter sampling used when `mode="time"`
+* `py=range(-dz/dhy,dz/dhy,64)` - ray parameter sampling used when `mode="time"`
+* `padz=0.0,padhx=0.0,padhy=0.0` - fractional padding in depth and offset to apply before applying the Fourier transform
+* `taperz=(0,0)` - beginning and end taper (fractional) in the z-direction (or t-direction) before transforming from `z-hx-hy` to `kz-khx-khy` or `t-hx-hy` to `f-khx-khy`
+* `taperhx=(0,0)` - beginning and end taper (fractional) in the hx-direction before transforming from `z-hx-hy` to `kz-khx-khy` or `t-hx-hy` to `f-khx-khy`
+* `taperhy=(0,0)` - beginning and end taper (fractional) in the hy-direction before transforming from `z-hx-hy` to `kz-khx-khy` or `t-hx-hy` to `f-khx-khy`
+* `taperkz=(0,0)` - beginning and end taper (fractional) in the kz-direction (or frequency) before sampling
+* `taperkhx=(0,0)` - beginning and end taper (fractional) in the khx-direction before sampling
+* `taperkhy=(0,0)` - beginning and end taper (fractional) in the khy-direction before sampling
+
+# Notes
+
+* If the mode is "time", then `padz` is the padding for the time dimension, `dz` is the time sampling interval, `taperz` is the taper for the time dimension and `tapkerkz` is the taper for frequency.
+* For mode="depth", typically `theta` needs to cover both positive and negative angles and must be < 90 degrees.
+* For mode="depth", the incidence angle depends on the geologic `dip` and `azimuth` (see 3-D Seismic Imaging by Prof. Biondo Biondi, Chapter 6). If not provided, this dependency is ignored.
+* Typically, geologic `dip` ∈ [0, 90] degrees and `azimuth` ∈ [0, 360] degrees, which is different from `theta` and `phi` ranges.
+* For mode="depth", only scalar `dip` and `azimuth` are supported, which means the same geologic angle is applied across the entire depth. This is a simplification and may not be accurate for complex geologies.
+"""
+function JopSlantStack3D(
+        dom::JetAbstractSpace{T};
+        theta = collect(-45.0:1.0:45.0),
+        phi = collect(0.0:45.0:135.0),
+        px = nothing,
+        py = nothing,
+        dz = 1.0,
+        dhx = 1.0,
+        dhy = 1.0,
+        hx0 = 0.0,
+        hy0 = 0.0,
+        padz = 0.0,
+        padhx = 0.0,
+        padhy = 0.0,
+        taperz = (0,0),
+        taperhx = (0,0),
+        taperhy = (0,0),
+        taperkz = (0,0),
+        taperkhx = (0,0),
+        taperkhy = (0,0),
+        dip = 0.0,
+        azimuth = 0.0,
+        mode = "depth") where {T}
+    mode ∈ ("depth", "time") || error("expected mode to be either 'depth' or 'time', got mode=$(mode)")
+    dz < 0.0 && error("expected dz>0.0, got dz=$(dz)")
+    dhx < 0.0 && error("expected dhx>0.0, got dhx=$(dhx)")
+    dhy < 0.0 && error("expected dhy>0.0, got dhy=$(dhy)")
+    mode == "depth" && any(abs.(theta) .>= 90.0) && error("for mode='depth', all theta values must be < 90 degrees in absolute value")
+
+    nz,nhx,nhy = size(dom)
+
+    # kz
+    nzfft = nextprod([2,3,5,7], round(Int, nz*(1 + padz)))
+    kz = rfftfreq(nzfft, 2*pi / dz)
+
+    # khx
+    nhxfft = nextprod([2,3,5,7], round(Int, nhx*(1+padhx)))
+    khx = fftfreq(nhxfft, 2*pi / dhx)
+
+    # khy
+    nhyfft = nextprod([2,3,5,7], round(Int, nhy*(1+padhy)))
+    khy = fftfreq(nhyfft, 2*pi / dhy)
+
+    # tan(theta), phi - used for mode=="depth"
+    tant = @. tan(deg2rad(theta))
+    phi = @. deg2rad(phi)
+
+    dip = deg2rad(dip)
+    azimuth = deg2rad(azimuth)
+
+    # p - used for both modes
+    if mode == "depth"
+        px = - tant
+        py = phi
+    elseif px === nothing || py === nothing
+        (px === nothing) && (px = collect(range(-dz/dhx, dz/dhx; length=64)))
+        (py === nothing) && (py = collect(range(-dz/dhy, dz/dhy; length=64)))
+    end
+
+    # conversions
+    kz,khx,khy,tant,phi,px,py = map(x->convert(Array{T,1}, x), (kz,khx,khy,tant,phi,px,py))
+    nzfft,nhxfft,nhyfft = map(x->convert(Int64, x), (nzfft,nhxfft,nhyfft))
+    hx0 = T(hx0)
+    hy0 = T(hy0)
+    dip = T(dip)
+    azimuth = T(azimuth)
+
+    # tapers
+    TX = JopTaper(dom, (1,2,3), (taperz[1],taperhx[1],taperhy[1]), (taperz[2],taperhx[2],taperhy[2]))
+    TK = JopTaper(JetSpace(Complex{eltype(dom)},div(nzfft,2)+1, nhxfft, nhyfft), (1,2,3), (taperkz[1], taperkhx[1], taperkhy[1]), (taperkz[2], taperkhx[2], taperkhy[2]), mode=(:normal,:fftshift,:fftshift))
+
+    JopLn(dom = dom, rng = JetSpace(T, nz, length(px), length(py)), df! = JopSlantStack3D_df!, df′! = JopSlantStack3D_df′!,
+        s = (;mode, nzfft, nhxfft, nhyfft, kz, khx, khy, px, py, hx0, hy0, dip, azimuth, TX, TK))
+end
+export JopSlantStack3D
+
+function JopSlantStack3D_df!(d::AbstractArray{T,3}, m::AbstractArray{T,3}; mode, nzfft, nhxfft, nhyfft, kz, khx, khy, px, py, hx0, hy0, dip, azimuth, TX, TK, kwargs...) where {T}
+    nhx, nhy, npx, npy, dhx, dhy = size(m,2), size(m,3), length(px), length(py), abs(khx[2]-khx[1]), abs(khy[2]-khy[1])
+    nz = size(d, 1)
+
+    mpad = zeros(T, nzfft, nhxfft, nhyfft)
+    mpad[1:nz,1:nhx,1:nhy] = TX*m
+    M = TK*rfft(mpad)
+
+    if mode == "time"
+        compute_kh = slantstack_khxy_from_frequency
+    elseif dip == 0
+        compute_kh = slantstack_khxy_from_kz
+    else
+        compute_kh = slantstack_khxy_from_kz_geologic
+    end
+
+    is_out_of_bounds = (ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1)->(ikhx_m1 < 1 || ikhx_p1 > nhxfft || ikhy_m1 < 1 || ikhy_p1 > nhyfft)
+
+    D = zeros(eltype(M), size(M,1), npx, npy)
+    @threads for ipx = 1:npx
+        for ipy = 1:npy
+            num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy])
+            denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2)
+            for ikz = 1:div(nzfft,2)+1
+                ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy = compute_kh(ikz, ipx, ipy, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom)
+
+                is_out_of_bounds(ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1) && continue
+
+                if (ikhx_m1 == ikhx_p1) && (ikhy_m1 == ikhy_p1)
+                    D[ikz,ipx,ipy] = M[ikz,ikhx_m1,ikhy_m1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0))
+                    continue
+                end
+
+                ax_m1 = (khx[ikhx_p1] - _khx)/dhx
+                ax_p1 = 1 - ax_m1
+                ay_m1 = (khy[ikhy_p1] - _khy)/dhy
+                ay_p1 = 1 - ay_m1
+
+                d_p1_p1 = M[ikz,ikhx_p1,ikhy_p1]*exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0))
+                d_p1_m1 = M[ikz,ikhx_p1,ikhy_m1]*exp(-im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0))
+                d_m1_p1 = M[ikz,ikhx_m1,ikhy_p1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0))
+                d_m1_m1 = M[ikz,ikhx_m1,ikhy_m1]*exp(-im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0))
+
+                D[ikz,ipx,ipy] = ax_m1*ay_m1*d_m1_m1 + ax_m1*ay_p1*d_m1_p1 + ax_p1*ay_m1*d_p1_m1 + ax_p1*ay_p1*d_p1_p1
+            end
+        end
+    end
+
+    _d = brfft(D, nzfft, 1)
+    d .= _d[1:nz,1:npx,1:npy] ./ nzfft
+end
+
+function JopSlantStack3D_df′!(m::AbstractArray{T,3}, d::AbstractArray{T,3}; mode, nzfft, nhxfft, nhyfft, kz, khx, khy, px, py, hx0, hy0, dip, azimuth, TX, TK, kwargs...) where {T}
+    nhx, nhy, npx, npy, dhx, dhy = size(m,2), size(m,3), length(px), length(py), abs(khx[2]-khx[1]), abs(khy[2]-khy[1])
+    nz = size(d, 1)
+
+    dpad = zeros(T, nzfft, npx, npy)
+    dpad[1:nz,:,:] = d
+
+    D = (rfft(dpad, 1) ./ nzfft)
+
+    if mode == "time"
+        compute_kh = slantstack_khxy_from_frequency
+    elseif dip == 0
+        compute_kh = slantstack_khxy_from_kz
+    else
+        compute_kh = slantstack_khxy_from_kz_geologic
+    end
+
+    is_out_of_bounds = (ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1)->(ikhx_m1 < 1 || ikhx_p1 > nhxfft || ikhy_m1 < 1 || ikhy_p1 > nhyfft)
+
+    M = zeros(Complex{T}, div(nzfft,2)+1, nhxfft, nhyfft)
+    for ipx = 1:npx
+        for ipy = 1:npy
+            num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy])
+            denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2)
+            for ikz = 1:div(nzfft,2)+1
+                ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy = compute_kh(ikz, ipx, ipy, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom)
+
+                is_out_of_bounds(ikhx_m1, ikhx_p1, ikhy_m1, ikhy_p1) && continue
+
+                if (ikhx_m1 == ikhx_p1) && (ikhy_m1 == ikhy_p1)
+                    M[ikz,ikhx_p1,ikhy_p1] += D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0))
+                    continue
+                end
+
+                ax_m1 = (khx[ikhx_p1] - _khx)/dhx
+                ax_p1 = 1 - ax_m1
+                ay_m1 = (khy[ikhy_p1] - _khy)/dhy
+                ay_p1 = 1 - ay_m1
+
+                m_p1_p1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_p1]*hy0))
+                m_p1_m1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_p1]*hx0 + khy[ikhy_m1]*hy0))
+                m_m1_p1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_p1]*hy0))
+                m_m1_m1 = D[ikz,ipx,ipy]*exp(im*(khx[ikhx_m1]*hx0 + khy[ikhy_m1]*hy0))
+
+                M[ikz,ikhx_p1,ikhy_p1] += ax_p1*ay_p1*m_p1_p1
+                M[ikz,ikhx_p1,ikhy_m1] += ax_p1*ay_m1*m_p1_m1
+                M[ikz,ikhx_m1,ikhy_p1] += ax_m1*ay_p1*m_m1_p1
+                M[ikz,ikhx_m1,ikhy_m1] += ax_m1*ay_m1*m_m1_m1
+            end
+        end
+    end
+
+    m .= TX * (brfft(TK * M, nzfft)[1:nz,1:nhx,1:nhy])
+end
+
+@inline function slantstack_khxy_from_kz(ikz::Int64, ipx::Int64, ipy::Int64, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom)
+    _khx = kz[ikz]*px[ipx]*cos(py[ipy])
+    _khy = kz[ikz]*px[ipx]*sin(py[ipy])
+
+    ikhx_m1 = floor(Int64, _khx/khx[2]) + 1
+    ikhx_p1 = ceil(Int64, _khx/khx[2]) + 1
+    ikhy_m1 = floor(Int64, _khy/khy[2]) + 1
+    ikhy_p1 = ceil(Int64, _khy/khy[2]) + 1
+
+    ikhx_m1 = ikhx_m1 < 1 ? nhxfft + ikhx_m1 : ikhx_m1
+    ikhx_p1 = ikhx_p1 < 1 ? nhxfft + ikhx_p1 : ikhx_p1
+    ikhy_m1 = ikhy_m1 < 1 ? nhyfft + ikhy_m1 : ikhy_m1
+    ikhy_p1 = ikhy_p1 < 1 ? nhyfft + ikhy_p1 : ikhy_p1
+
+    ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy
+end
+
+@inline function slantstack_khxy_from_kz_geologic(ikz::Int64, ipx::Int64, ipy::Int64, px, py, kz, khx, khy, nhxfft, nhyfft, num, denom)
+    _khx = kz[ikz]*px[ipx]*(cos(py[ipy])/denom + sin(py[ipy]*num/denom))
+    _khy = kz[ikz]*px[ipx]*(sin(py[ipy])/denom - cos(py[ipy]*num/denom))
+
+    ikhx_m1 = floor(Int64, _khx/khx[2]) + 1
+    ikhx_p1 = ceil(Int64, _khx/khx[2]) + 1
+    ikhy_m1 = floor(Int64, _khy/khy[2]) + 1
+    ikhy_p1 = ceil(Int64, _khy/khy[2]) + 1
+
+    ikhx_m1 = ikhx_m1 < 1 ? nhxfft + ikhx_m1 : ikhx_m1
+    ikhx_p1 = ikhx_p1 < 1 ? nhxfft + ikhx_p1 : ikhx_p1
+    ikhy_m1 = ikhy_m1 < 1 ? nhyfft + ikhy_m1 : ikhy_m1
+    ikhy_p1 = ikhy_p1 < 1 ? nhyfft + ikhy_p1 : ikhy_p1
+
+    ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy
+end
+
+@inline function slantstack_khxy_from_frequency(iω::Int64, ipx::Int64, ipy::Int64, px, py, ω, khx, khy, nhxfft, nhyfft, num, denom)
+    _khx = px[ipx]*ω[iω]
+    _khy = py[ipy]*ω[iω]
+
+    ikhx_m1 = floor(Int64, _khx/khx[2]) + 1
+    ikhx_p1 = ceil(Int64, _khx/khx[2]) + 1
+    ikhy_m1 = floor(Int64, _khy/khy[2]) + 1
+    ikhy_p1 = ceil(Int64, _khy/khy[2]) + 1
+
+    ikhx_m1 = ikhx_m1 < 1 ? nhxfft + ikhx_m1 : ikhx_m1
+    ikhx_p1 = ikhx_p1 < 1 ? nhxfft + ikhx_p1 : ikhx_p1
+    ikhy_m1 = ikhy_m1 < 1 ? nhyfft + ikhy_m1 : ikhy_m1
+    ikhy_p1 = ikhy_p1 < 1 ? nhyfft + ikhy_p1 : ikhy_p1
+
+    ikhx_m1, ikhx_p1, _khx, ikhy_m1, ikhy_p1, _khy
+end
+
+"""
+    A = JopSlantStackShiftSum(dom[; dz=1.0, ...])
+
+where `A` is the 2D slant-stack operator mapping for `z-h` to `z-θ` (depth mode) or `t-h` to `tau-p` (time mode).
+The slant stacking is performed directly on the input domain by shifting and summing along the offset axis.
+The domain of the operator is `nz` x `nh` with precision T, `dz` is the depth spacing (or time interval),
+`h`, if provided, is the array of offsets (can be irregular). If not provided, a regular grid is assumed with same sampling as dz.
+The additional named optional arguments along with their default values are,
+
+* `mode="depth` - choose between "depth" and "time" to specify if the input domain is `z-h` or `t-h`
+* `theta=-45:1.0:45` - range of incidence angles used when `mode="depth"`
+* `p=range(-dz/max(diff(h)),dz/max(diff(h)),128)` - ray parameter sampling used when `mode="time"`
+* `taperz=(0,0)` - beginning and end taper (fractional) in the z-direction (or t-direction) before shifting and summing
+
+# Notes
+
+* If the mode is "time", then `dz` is the time sampling interval, `taperz` is the taper for the time dimension.
+* For mode="depth", typically theta needs to cover both positive and negative angles and must be < 90 degrees.
+"""
+function JopSlantStackShiftSum(
+        dom::JetAbstractSpace{T};
+        theta = collect(-45.0:1.0:45.0),
+        p = nothing,
+        dz = 1.0,
+        h = nothing,
+        taperz = (0,0),
+        mode = "depth") where {T}
+    mode ∈ ("depth", "time") || error("expected mode to be either 'depth' or 'time', got mode=$(mode)")
+    dz < 0.0 && error("expected dz>0.0, got dz=$(dz)")
+    h === nothing && (h = collect(0.0:dz:(size(dom,2)-1)*dz))
+    length(h) != size(dom,2) && error("length of h=$(length(h)) must match the number of offsets in the domain=$(size(dom,2))")
+    dhmax = maximum(diff(sort(h)))
+    mode == "depth" && any(abs.(theta) .>= 90.0) && error("for mode='depth', all theta values must be < 90 degrees in absolute value")
+
+    nz,nh = size(dom)
+
+    # tan(theta) - used for mode=="depth"
+    tant = @. tan(deg2rad(theta))
+
+    # p - used for both modes
+    if mode == "depth"
+        p = - tant
+    elseif p === nothing
+        p = collect(range(-dz/dhmax, dz/dhmax; length=128))
+    end
+
+    # conversions
+    h,p = map(x->convert(Array{T,1}, x), (h,p))
+    dz = T(dz)
+
+    # taper
+    TZ = JopTaper(dom, (1,), (taperz[1],), (taperz[2],))
+
+    JopLn(dom = dom, rng = JetSpace(T, nz, length(p)), df! = JopSlantStackShiftSum_df!, df′! = JopSlantStackShiftSum_df′!,
+        s = (;dz, h, p, TZ))
+end
+export JopSlantStackShiftSum
+
+function JopSlantStackShiftSum_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; dz, h, p, TZ, kwargs...) where {T}
+    nz, nh, np = size(m,1), size(m,2), length(p)
+
+    d .= 0
+    mtap = TZ * m
+    @threads for ip = 1:np
+        holder = zeros(T, nz)
+        for ih = 1:nh
+            shift = + h[ih] * p[ip] / dz
+            if abs(shift) < nz
+                holder=_shift_forward(@view(mtap[:,ih]), shift)
+                d[:,ip] .+= holder
+            end
+        end
+    end
+    d
+end
+
+function JopSlantStackShiftSum_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; dz, h, p, TZ, kwargs...) where {T}
+    nz, nh, np = size(m,1), size(m,2), length(p)
+
+    mtap = zeros(T, nz, nh)
+    @threads for ih = 1:nh
+        holder = zeros(T, nz)
+        for ip = 1:np
+            shift = + h[ih] * p[ip] / dz
+            if abs(shift) < nz
+                holder=_shift_adjoint(@view(d[:,ip]), shift)
+                @views mtap[:,ih] .+= holder
+            end
+        end
+    end
+    m = TZ' * mtap
+    m
+end
+
+
+
+
+"""
+    A = JopSlantStackShiftSum3D(dom[; dz=1.0, ...])
+
+where `A` is the 3D slant-stack operator mapping for `z-hx-hy` to `z-θ-ϕ` (depth mode) or `t-hx-hy` to `tau-px-py` (time mode).
+The slant stacking is performed directly on the input domain by shifting and summing along the offset axis.
+The domain of the operator is `nz` x `nh` (for irregular) or `nz` x `nhx` x `nhy` (for regular) with precision T, `dz` is the depth spacing (or time interval),
+`hx` and `hy`, if provided, are the arrays of offsets (can be irregular). If not provided, a regular grid is assumed with same sampling as dz.
+The additional named optional arguments along with their default values are,
+
+* `mode="depth` - choose between "depth" and "time" to specify if the input domain is `z-hx-hy` or `t-hx-hy`
+* `theta=-45:1.0:45` - range of incidence angles used when `mode="depth"`
+* `phi=0:45.0:135` - range of image azimuth angles used when `mode="depth"`
+* `dip=0.0` - geologic dip used when `mode="depth"`
+* `azimuth=0.0` - geologic azimuth used when `mode="depth"`
+* `px=range(-dz/max(diff(hx)),dz/max(diff(hx)),64)` - ray parameter sampling used when `mode="time"`
+* `py=range(-dz/max(diff(hy)),dz/max(diff(hy)),64)` - ray parameter sampling used when `mode="time"`
+* `taperz=(0,0)` - beginning and end taper (fractional) in the z-direction (or t-direction) before shifting and summing
+
+# Notes
+
+* If the mode is "time", then `dz` is the time sampling interval, `taperz` is the taper for the time dimension.
+* For mode="depth", typically `theta` needs to cover both positive and negative angles and must be < 90 degrees.
+* For mode="depth", the incidence angle depends on the geologic `dip` and `azimuth` (see 3-D Seismic Imaging by Prof. Biondo Biondi, Chapter 6). If not provided, this dependency is ignored.
+* Typically, geologic `dip` ∈ [0, 90] degrees and `azimuth` ∈ [0, 360] degrees, which is different from `theta` and `phi` ranges.
+"""
+function JopSlantStackShiftSum3D(
+        dom::JetAbstractSpace{T};
+        theta = collect(-45.0:1.0:45.0),
+        phi = collect(0.0:45.0:135.0),
+        px = nothing,
+        py = nothing,
+        dz = 1.0,
+        hx = nothing,
+        hy = nothing,
+        taperz = (0,0),
+        dip = 0.0,
+        azimuth = 0.0,
+        mode = "depth") where {T}
+    mode ∈ ("depth", "time") || error("expected mode to be either 'depth' or 'time', got mode=$(mode)")
+    dz < 0.0 && error("expected dz>0.0, got dz=$(dz)")
+    if mode == "time" || all(iszero, dip)
+        dip, azimuth = 0, 0
+    end
+    (typeof(dip) != typeof(azimuth)) && error("dip and azimuth should be of the same type")
+    isa(dip, Array) && (length(dip) != length(azimuth)) && error("if dip is an array, its length=$(length(dip)) should match the length of azimuth=$(length(azimuth))")
+    isa(dip, Array) && (length(dip) != size(dom, 1)) && error("if dip is an array, its length=$(length(dip)) should match the first dimension of the domain=$(size(dom, 1))")
+
+    nd = ndims(dom) # if nd = 2, the offsets will be considered as totally irregular
+    nd ∉ (2,3) && error("expected domain to be either 2D or 3D, got ndims=$(nd)")
+    (nd == 2) && (hx === nothing || hy === nothing) && error("for 2D irregular (non-lattice) offsets, hx and hy should be provided") 
+    ndy = (nd == 2 ? 2 : 3)
+    hx === nothing && (hx = collect(0.0:dz:(size(dom,2)-1)*dz))
+    hy === nothing && (hy = collect(0.0:dz:(size(dom,3)-1)*dz))
+    length(hx) != size(dom,2) && error("length of hx=$(length(hx)) must match the number of x-offsets in the domain=$(size(dom,2))")
+    length(hy) != size(dom,ndy) && error("length of hy=$(length(hy)) must match the number of y-offsets in the domain=$(size(dom,ndy))")
+    dhxmax = length(hx) > 1 ? maximum(diff(sort(hx))) : 1.0
+    dhymax = length(hy) > 1 ? maximum(diff(sort(hy))) : 1.0
+    mode == "depth" && any(abs.(theta) .>= 90.0) && error("for mode='depth', theta values must be < 90 degrees in absolute value")
+    mode == "depth" && (any(phi .> 180.0) || any(phi .< 0.0)) && error("for mode='depth', phi values must be between 0 and 180 degrees")
+
+    nz = size(dom, 1)
+    nhx = size(dom, 2)
+    nhy = size(dom, ndy)
+
+    # tan(theta), phi - used for mode=="depth"
+    tant = @. tan(deg2rad(theta))
+    phi = @. deg2rad(phi)
+
+    dip = @. deg2rad(dip)
+    azimuth = @. deg2rad(azimuth)
+
+    # p - used for both modes
+    if mode == "depth"
+        px = - tant
+        py = phi
+    elseif px === nothing || py === nothing
+        (px === nothing) && (px = collect(range(-dz/dhxmax, dz/dhxmax; length=64)))
+        (py === nothing) && (py = collect(range(-dz/dhymax, dz/dhymax; length=64)))
+    end
+
+    # conversions
+    hx,hy,px,py = map(x->convert(Array{T,1}, x), (hx,hy,px,py))
+    if !isa(dip, Array)
+        dip = T(dip)
+        azimuth = T(azimuth)
+    else
+        dip, azimuth = map(x->convert(Array{T,1}, x), (dip, azimuth))
+    end
+    dz = T(dz)
+
+    # taper
+    TZ = JopTaper(dom, (1,), (taperz[1],), (taperz[2],))
+
+    JopLn(dom = dom, rng = JetSpace(T, nz, length(px), length(py)), df! = JopSlantStackShiftSum3D_df!, df′! = JopSlantStackShiftSum3D_df′!,
+        s = (;mode, dz, hx, hy, px, py, dip, azimuth, TZ))
+end
+export JopSlantStackShiftSum3D
+
+function JopSlantStackShiftSum3D_df!(d::AbstractArray{T,3}, m::AbstractArray{T,2}; dip, azimuth, mode, dz, hx, hy, px, py, TZ, kwargs...) where {T}
+    if isa(dip, Array)
+        JopSlantStackShiftSum3D_df_vector!(d, m; dip, azimuth, mode, dz, hx, hy, px, py, TZ)
+    else
+        JopSlantStackShiftSum3D_df_scalar!(d, m; dip, azimuth, mode, dz, hx, hy, px, py, TZ)
+    end
+end
+
+function JopSlantStackShiftSum3D_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,3}; dip, azimuth, mode, dz, hx, hy, px, py, TZ, kwargs...) where {T}
+    if isa(dip, Array)
+        JopSlantStackShiftSum3D_df_vector′!(m, d; dip, azimuth, mode, dz, hx, hy, px, py, TZ, kwargs...)
+    else
+        JopSlantStackShiftSum3D_df_scalar′!(m, d; dip, azimuth, mode, dz, hx, hy, px, py, TZ, kwargs...)
+    end
+end
+
+function JopSlantStackShiftSum3D_df!(d::AbstractArray{T,3}, m::AbstractArray{T,3}; dip, azimuth, mode, dz, hx, hy, px, py, TZ, kwargs...) where {T}
+    if isa(dip, Array)
+        JopSlantStackShiftSum3D_df_vector!(d, m; dip, azimuth, mode, dz, hx, hy, px, py, TZ)
+    else
+        JopSlantStackShiftSum3D_df_scalar!(d, m; dip, azimuth, mode, dz, hx, hy, px, py, TZ)
+    end
+end
+
+function JopSlantStackShiftSum3D_df′!(m::AbstractArray{T,3}, d::AbstractArray{T,3}; dip, azimuth, mode, dz, hx, hy, px, py, TZ, kwargs...) where {T}
+    if isa(dip, Array)
+        JopSlantStackShiftSum3D_df_vector′!(m, d; dip, azimuth, mode, dz, hx, hy, px, py, TZ, kwargs...)
+    else
+        JopSlantStackShiftSum3D_df_scalar′!(m, d; dip, azimuth, mode, dz, hx, hy, px, py, TZ, kwargs...)
+    end
+end
+
+function JopSlantStackShiftSum3D_df_scalar!(d::AbstractArray{T,3}, m::AbstractArray{T,2}; dip::Real, azimuth::Real, mode, dz, hx, hy, px, py, TZ) where {T}
+    nz, nh, npx, npy = size(m,1), size(m,2), length(px), length(py)
+
+    if mode == "time"
+        compute_shift = slantstack_shift_px_py
+    elseif dip == 0
+        compute_shift = slantstack_shift_theta_phi
+    else
+        compute_shift = slantstack_shift_theta_phi_geologic
+    end
+
+    d .= 0
+    mtap = TZ * m
+    @threads for ipx = 1:npx
+        holder = zeros(T, nz)
+        for ipy = 1:npy
+            num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy])
+            denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2)
+            for ih = 1:nh
+                shift = compute_shift(1, ih, ih, ipx, ipy, hx, hy, px, py, dz, num, denom)
+                if abs(shift) < nz
+                    holder=_shift_forward(@view(mtap[:,ih]), shift)
+                    d[:,ipx,ipy] .+= holder
+                end
+            end
+        end
+    end
+    d
+end
+
+function JopSlantStackShiftSum3D_df_scalar!(d::AbstractArray{T,3}, m::AbstractArray{T,3}; dip::Real, azimuth::Real, mode, dz, hx, hy, px, py, TZ) where {T}
+    nz, nhx, nhy, npx, npy = size(m,1), size(m,2), size(m,3), length(px), length(py)
+
+    if mode == "time"
+        compute_shift = slantstack_shift_px_py
+    elseif dip == 0
+        compute_shift = slantstack_shift_theta_phi
+    else
+        compute_shift = slantstack_shift_theta_phi_geologic
+    end
+
+    d .= 0
+    mtap = TZ * m
+    @threads for ipx = 1:npx
+        holder = zeros(T, nz)
+        for ipy = 1:npy
+            num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy])
+            denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2)
+            for ihx = 1:nhx
+                for ihy = 1:nhy
+                    shift = compute_shift(1, ihx, ihy, ipx, ipy, hx, hy, px, py, dz, num, denom)
+                    if abs(shift) < nz
+                        holder=_shift_forward(@view(mtap[:,ihx,ihy]), shift)
+                        d[:,ipx,ipy] .+= holder
+                    end
+                end
+            end
+        end
+    end
+    d
+end
+
+function JopSlantStackShiftSum3D_df_vector!(d::AbstractArray{T,3}, m::AbstractArray{T,2}; dip::AbstractArray{T,1}, azimuth::AbstractArray{T,1}, mode, dz, hx, hy, px, py, TZ) where {T}
+    nz, nh, npx, npy = size(m,1), size(m,2), length(px), length(py)
+
+    compute_shift = slantstack_shift_theta_phi_geologic
+    
+    d .= 0
+    mtap = TZ * m
+    @threads for ipx = 1:npx
+        holder = zeros(T, nz)
+        for ipy = 1:npy
+            num = tan.(dip).^2 .* sin.(azimuth .- py[ipy]) .* cos.(azimuth .- py[ipy])
+            denom = sqrt.(1 .+ (tan.(dip) .* sin.(azimuth .- py[ipy])).^2)
+            for ih = 1:nh
+                for iz = 1:nz
+                    shift = compute_shift(iz, ih, ih, ipx, ipy, hx, hy, px, py, dz, num, denom)
+                    if abs(shift) < nz 
+                        holder=_shift_forward(@view(mtap[:,ih]), shift)
+                        d[iz,ipx,ipy] += holder[iz]
+                    end
+                end
+            end
+        end
+    end
+    d
+end
+
+function JopSlantStackShiftSum3D_df_vector!(d::AbstractArray{T,3}, m::AbstractArray{T,3}; dip::AbstractArray{T,1}, azimuth::AbstractArray{T,1}, mode, dz, hx, hy, px, py, TZ) where {T}
+    nz, nhx, nhy, npx, npy = size(m,1), size(m,2), size(m,3), length(px), length(py)
+
+    compute_shift = slantstack_shift_theta_phi_geologic
+    
+    d .= 0
+    mtap = TZ * m
+    @threads for ipx = 1:npx
+        holder = zeros(T, nz)
+        for ipy = 1:npy
+            num = tan.(dip).^2 .* sin.(azimuth .- py[ipy]) .* cos.(azimuth .- py[ipy])
+            denom = sqrt.(1 .+ (tan.(dip) .* sin.(azimuth .- py[ipy])).^2)
+            for ihx = 1:nhx
+                for ihy = 1:nhy
+                    for iz = 1:nz
+                        shift = compute_shift(iz, ihx, ihy, ipx, ipy, hx, hy, px, py, dz, num, denom)
+                        if abs(shift) < nz 
+                            holder=_shift_forward(@view(mtap[:,ihx,ihy]), shift)
+                            d[iz,ipx,ipy] += holder[iz]
+                        end
+                    end
+                end
+            end
+        end
+    end
+    d
+end
+
+function JopSlantStackShiftSum3D_df_scalar′!(m::AbstractArray{T,2}, d::AbstractArray{T,3}; dip::Real, azimuth::Real, mode, dz, hx, hy, px, py, TZ, kwargs...) where {T}
+    nz, nh, npx, npy = size(m,1), size(m,2), length(px), length(py)
+
+    if mode == "time"
+        compute_shift = slantstack_shift_px_py
+    elseif dip == 0
+        compute_shift = slantstack_shift_theta_phi
+    else
+        compute_shift = slantstack_shift_theta_phi_geologic
+    end
+
+    mtap = zeros(T, nz, nh)
+    @threads for ih = 1:nh
+        acc = zeros(T, nz)
+        holder = zeros(T, nz)
+        for ipx = 1:npx
+            for ipy = 1:npy
+                num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy])
+                denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2)
+                shift = compute_shift(1, ih, ih, ipx, ipy, hx, hy, px, py, dz, num, denom)
+                if abs(shift) < nz
+                    holder=_shift_adjoint(@view(d[:,ipx,ipy]), shift)
+                    acc .+= holder
+                end
+            end
+        end
+        mtap[:,ih] .= acc
+    end
+    m = TZ' * mtap
+    m
+end
+
+function JopSlantStackShiftSum3D_df_scalar′!(m::AbstractArray{T,3}, d::AbstractArray{T,3}; dip::Real, azimuth::Real, mode, dz, hx, hy, px, py, TZ, kwargs...) where {T}
+    nz, nhx, nhy, npx, npy = size(m,1), size(m,2), size(m,3), length(px), length(py)
+
+    if mode == "time"
+        compute_shift = slantstack_shift_px_py
+    elseif dip == 0
+        compute_shift = slantstack_shift_theta_phi
+    else
+        compute_shift = slantstack_shift_theta_phi_geologic
+    end
+
+    mtap = zeros(T, nz, nhx, nhy)
+    @threads for ihx = 1:nhx
+        for ihy = 1:nhy
+            acc = zeros(T, nz)
+            holder = zeros(T, nz)
+            for ipx = 1:npx
+                for ipy = 1:npy
+                    num = tan(dip)^2 *sin(azimuth - py[ipy])*cos(azimuth - py[ipy])
+                    denom = sqrt(1 + (tan(dip)*sin(azimuth - py[ipy]))^2)
+                    shift = compute_shift(1, ihx, ihy, ipx, ipy, hx, hy, px, py, dz, num, denom)
+                    if abs(shift) < nz
+                        holder=_shift_adjoint(@view(d[:,ipx,ipy]), shift)
+                        acc .+= holder
+                    end
+                end
+            end
+            mtap[:,ihx,ihy] .= acc
+        end
+    end
+    m = TZ' * mtap
+    m
+end
+
+function JopSlantStackShiftSum3D_df_vector′!(m::AbstractArray{T,2}, d::AbstractArray{T,3}; dip::AbstractArray{T,1}, azimuth::AbstractArray{T,1}, mode, dz, hx, hy, px, py, TZ, kwargs...) where {T}
+    nz, nh, npx, npy = size(m,1), size(m,2), length(px), length(py)
+
+    compute_shift = slantstack_shift_theta_phi_geologic
+
+    mtap = zeros(T, nz, nh)
+    @threads for ih = 1:nh
+        acc = zeros(T, nz)
+        holder = zeros(T, nz)
+        for ipx = 1:npx
+            for ipy = 1:npy
+                num = tan.(dip).^2 .* sin.(azimuth .- py[ipy]) .* cos.(azimuth .- py[ipy])
+                denom = sqrt.(1 .+ (tan.(dip) .* sin.(azimuth .- py[ipy])).^2)
+                for iz = 1:nz
+                    shift = compute_shift(iz, ih, ih, ipx, ipy, hx, hy, px, py, dz, num, denom)
+                    if abs(shift) < nz
+                        fill!(holder, zero(T))
+                        holder[iz] = d[iz,ipx,ipy]
+                        acc=_shift_adjoint(holder, shift)
+                        mtap[:,ih] .+= acc
+                    end
+                end
+            end
+        end
+    end
+    m = TZ' * mtap
+    m
+end
+
+function JopSlantStackShiftSum3D_df_vector′!(m::AbstractArray{T,3}, d::AbstractArray{T,3}; dip::AbstractArray{T,1}, azimuth::AbstractArray{T,1}, mode, dz, hx, hy, px, py, TZ, kwargs...) where {T}
+    nz, nhx, nhy, npx, npy = size(m,1), size(m,2), size(m,3), length(px), length(py)
+
+    compute_shift = slantstack_shift_theta_phi_geologic
+
+    mtap = zeros(T, nz, nhx, nhy)
+    @threads for ihx = 1:nhx
+        for ihy = 1:nhy
+            acc = zeros(T, nz)
+            holder = zeros(T, nz)
+            for ipx = 1:npx
+                for ipy = 1:npy
+                    num = tan.(dip).^2 .* sin.(azimuth .- py[ipy]) .* cos.(azimuth .- py[ipy])
+                    denom = sqrt.(1 .+ (tan.(dip) .* sin.(azimuth .- py[ipy])).^2)
+                    for iz = 1:nz
+                        shift = compute_shift(iz, ihx, ihy, ipx, ipy, hx, hy, px, py, dz, num, denom)
+                        if abs(shift) < nz
+                            fill!(holder, zero(T))
+                            holder[iz] = d[iz,ipx,ipy]
+                            acc=_shift_adjoint(holder, shift)
+                            mtap[:,ihx,ihy] .+= acc
+                        end
+                    end
+                end
+            end
+        end
+    end
+    m = TZ' * mtap
+    m
+end
+
+@inline function slantstack_shift_px_py(iz::Int64, ihx::Int64, ihy::Int64, ipx::Int64, ipy::Int64, hx, hy, px, py, dz, num::Real, denom::Real)
+    # for time domain, the shift is given by (hx*px + hy*py) / dz
+    shift = + (hx[ihx] * px[ipx] + hy[ihy] * py[ipy]) / dz
+    shift
+end
+
+@inline function slantstack_shift_theta_phi(iz::Int64, ihx::Int64, ihy::Int64, ipx::Int64, ipy::Int64, hx, hy, px, py, dz, num::Real, denom::Real)
+    # for depth domain ignoring geologic dip, the shift is given by (-hx*tan(θ)*cos(ϕ) - hy*tan(θ)*sin(ϕ)) / dz
+    shift = + (hx[ihx] * px[ipx] * cos(py[ipy]) + hy[ihy] * px[ipx] * sin(py[ipy])) / dz
+    shift
+end
+
+@inline function slantstack_shift_theta_phi_geologic(iz::Int64, ihx::Int64, ihy::Int64, ipx::Int64, ipy::Int64, hx, hy, px, py, dz, num::Real, denom::Real)
+    # for depth domain accounting for geologic dip, the shift is given by (-hx*tan(θ)*(cos(ϕ)/denom + sin(ϕ)*num/denom) - hy*tan(θ)*(sin(ϕ)/denom - cos(ϕ)*num/denom)) / dz
+    shift = + hx[ihx] * px[ipx] * (cos(py[ipy]) / denom + sin(py[ipy]) * num / denom) + hy[ihy] * px[ipx] * (sin(py[ipy]) / denom - cos(py[ipy]) * num / denom)
+    shift / dz
+end
+
+@inline function slantstack_shift_theta_phi_geologic(iz::Int64, ihx::Int64, ihy::Int64, ipx::Int64, ipy::Int64, hx, hy, px, py, dz, num::AbstractArray{T,1}, denom::AbstractArray{T,1}) where {T}
+    # for depth domain accounting for geologic dip, the shift is given by (-hx*tan(θ)*(cos(ϕ)/denom + sin(ϕ)*num/denom) - hy*tan(θ)*(sin(ϕ)/denom - cos(ϕ)*num/denom)) / dz
+    shift = + hx[ihx] * px[ipx] * (cos(py[ipy]) / denom[iz] + sin(py[ipy]) * num[iz] / denom[iz]) + hy[ihy] * px[ipx] * (sin(py[ipy]) / denom[iz] - cos(py[ipy]) * num[iz] / denom[iz])
+    shift / dz
+end
+
+# helper functions to build a compact sinc kernel for forward and adjoint shifting in 1D
+function _sinc_kernel(δ, N)
+    n = -(N÷2):(N÷2)
+    h = sinc.(n .- δ)
+    h .* hann(length(h))
+end
+
+function _shift_forward(x, shift, N = 7)
+    ishift = round(Int, shift)
+    frac = shift - ishift
+    h = _sinc_kernel(frac, N)
+    L = length(h) ÷ 2
+    x_int = circshift(x, ishift)
+    x_frac = conv(x_int, h)
+    x_frac[L+1 : L+length(x)]
+end
+
+function _shift_adjoint(x, shift, N = 7)
+    ishift = round(Int, shift)
+    frac = shift - ishift
+    h = reverse(_sinc_kernel(frac, N))
+    n = length(x)
+    m = length(h)
+    L = m ÷ 2
+    x_pad = zeros(eltype(x), n+m-1)
+    x_pad[L+1 : L+n] .= x
+    x_frac = conv(x_pad, h)
+    circshift(@view(x_frac[m : m+n-1]), -ishift)
+end
\ No newline at end of file
diff --git a/test/jop_slantstack.jl b/test/jop_slantstack.jl
index 83d4c2a..e88ae39 100644
--- a/test/jop_slantstack.jl
+++ b/test/jop_slantstack.jl
@@ -4,9 +4,6 @@ using JetPackTransforms, Jets, Test, LinearAlgebra
 @testset "JetSlantStack, dot product test" for T in (Float32, Float64)
     A = JopSlantStack(JetSpace(T, 64, 128); dz=10.0, dh=10.0, h0=-1000.0)
 
-    m = rand(domain(A))
-    d = A*m
-
     lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
     @test isapprox(lhs,rhs,rtol=1e-4)
 
@@ -29,9 +26,6 @@ end
 @testset "JetSlantStack, dot product test" for T in (Float32, Float64)
     A = JopSlantStack(JetSpace(T, 64, 128); dz=0.004, dh=10.0, h0=-1000.0, mode="time")
 
-    m = rand(domain(A))
-    d = A*m
-
     lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
     @test isapprox(lhs,rhs,rtol=1e-4)
 
@@ -51,9 +45,8 @@ end
 
 # depth-time parity test
 @testset "JetSlantStack, parity" begin
-    theta = collect(-45:0.5:45)
+    theta = collect(-45:1.0:45)
     # the ray parameters that would give the same results as theta is
-    # p = -tan(θ)
     p = @. -tan(deg2rad(theta))
     
     Ad = JopSlantStack(JetSpace(Float64, 128, 129); mode="depth", theta=theta, p=p, dz=0.01, dh=0.01, h0=-0.64, padz=1, padh=1, taperkz = (0.5,0.5), taperkh = (0.5,0.5))
@@ -67,12 +60,246 @@ end
     dd = Ad*m
     dt = At*m
 
-    md = Ad'*dd
-    mt = At'*dt
+    err_d = maximum(abs.(dd .- dt)) / maximum(abs.(dd))
+    @show err_d
+    @test err_d < 1e-7
+end
+
+@testset "JetSlantStackShiftSum, dot product test" for T in (Float32, Float64), mode in ("depth", "time")
+    A = JopSlantStackShiftSum(JetSpace(T, 64, 128); dz=10.0, mode=mode, h = 10.0 .* rand(128) .- 5, taperz=(0.3,0.3))
+    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
+    @test isapprox(lhs,rhs,rtol=1e-4)
+end
+
+@testset "JetSlantStackShiftSum, parity" begin
+    theta = collect(-45:1.0:45)
+    # the ray parameters that would give the same results as theta is
+    p = @. -tan(deg2rad(theta))
+    
+    Ad = JopSlantStackShiftSum(JetSpace(Float64, 128, 129); mode="depth", theta=theta, p=p, dz=0.01, h=collect(-0.64:0.01:0.64))
+    At = JopSlantStackShiftSum(JetSpace(Float64, 128, 129); mode="time" , theta=theta, p=p, dz=0.01, h=collect(-0.64:0.01:0.64))
+    m = zeros(domain(Ad))
+    for i = 1:129
+        m[div(i,4)+20,i] = 1.0
+    end
+    m[16,64] = 1
+    m[64,96] = 1
+    dd = Ad*m
+    dt = At*m
 
     err_d = maximum(abs.(dd .- dt)) / maximum(abs.(dd))
-    err_m = maximum(abs.(md .- mt)) / maximum(abs.(md))
-    @show err_d, err_m
+
+    @show err_d
     @test err_d < 1e-7
-    @test err_m < 1e-7
 end
+
+@testset "JetSlantStack vs JetSlantStackShiftSum, parity" begin
+    theta = collect(-45:1.0:45)
+    A1 = JopSlantStack(JetSpace(Float64, 128, 129); mode="depth", theta=theta, dz=0.01, h0=-0.64, dh=0.01)
+    A2 = JopSlantStackShiftSum(JetSpace(Float64, 128, 129); mode="depth" , theta=theta, dz=0.01, h=collect(-0.64:0.01:0.64))
+    m = zeros(domain(A1))
+    for i = 1:129
+        m[div(i,4)+20,i] = 1.0
+    end
+    m[16,64] = 1
+    m[64,96] = 1
+    d1 = A1*m
+    d2 = A2*m
+
+    # check cosine similarity of the results since the two operators are not exactly the same
+    similarity = dot(d1, d2) / (norm(d1) * norm(d2))
+
+    @show similarity
+    @test similarity > 0.95
+end
+
+@testset "JetSlantStack3D, dot product test, T = $(T), mode = $(mode), dip = $(dip)" for T in (Float32, Float64), mode in ("depth", "time"), (dip, azimuth) in ((0.0,0.0), (30.0, 60.0))
+    A = JopSlantStack3D(JetSpace(T, 32, 8, 5); dz=10.0, dhx=10.0, dhy=5.0, hx0=-1000.0, hy0=-500.0, mode=mode, dip=dip, azimuth=azimuth, taperz=(0.3,0.3), taperhx=(0.3,0.3), taperhy=(0.3,0.3), taperkz=(0.3,0.3), taperkhx=(0.3,0.3), taperkhy=(0.3,0.3))
+
+    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
+    @test isapprox(lhs,rhs,rtol=1e-4)
+end
+
+@testset "JetSlantStack3D, time-depth parity" begin
+    theta = collect(-45:1.0:45)
+    phi = [0.0, 90.0]
+    
+    # the ray parameters that would give the same results as theta/phi are
+    px = @. -tan(deg2rad(theta))
+    py = [0.0]
+    
+    Ad = JopSlantStack3D(JetSpace(Float64, 64, 8, 8); mode="depth", theta=theta, phi = phi, px=px, py=py, dz=0.01, dhx=0.01, dhy=0.01, hx0=-0.64, hy0=-0.64)
+    At1 = JopSlantStack3D(JetSpace(Float64, 64, 8, 8); mode="time" , theta=theta, phi = phi, px=px, py=py, dz=0.01, dhx=0.01, dhy=0.01, hx0=-0.64, hy0=-0.64)
+    At2 = JopSlantStack3D(JetSpace(Float64, 64, 8, 8); mode="time" , theta=theta, phi = phi, px=py, py=px, dz=0.01, dhx=0.01, dhy=0.01, hx0=-0.64, hy0=-0.64)
+    m = zeros(domain(Ad))
+    for i = 1:8
+        m[div(i,4)+20,i,:] .= 1.0
+        m[div(i,2)+20,:,i] .= 0.5
+    end
+    m[16,4,4] = 1
+    m[32,5,7] = 1
+    dd = Ad*m
+    dt1 = At1*m
+    dt2 = At2*m
+
+    err_d1 = maximum(abs.(dd[:,:,1] .- dt1[:,:,1])) / maximum(abs.(dd[:,:,1]))
+    err_d2 = maximum(abs.(dd[:,:,2] .- dt2[:,1,:])) / maximum(abs.(dd[:,:,2]))
+
+    @show err_d1, err_d2
+    @test err_d1 < 1e-7
+    @test err_d2 < 1e-7
+end
+
+@testset "JetSlantStackShiftSum3D, dot product test, T = $(T), mode = $(mode), dip = $(dip)" for T in (Float32, Float64), mode in ("depth", "time"), (dip, azimuth) in ((0.0,0.0), (30.0, 60.0), ([0.0], [0.0]))
+    if isa(dip, Array)
+        dip = 90 .* rand(32)
+        azimuth = 360 .* rand(32)
+    end
+    A = JopSlantStackShiftSum3D(JetSpace(T, 32, 64); dz=10.0, mode=mode, hx = 10.0 .* rand(64) .- 5, hy = 10.0 .* rand(64) .- 5,  dip=dip, azimuth=azimuth, taperz=(0.3,0.3))
+    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
+    @test isapprox(lhs,rhs,rtol=1e-4)
+
+    A = JopSlantStackShiftSum3D(JetSpace(T, 32, 8, 8); dz=10.0, mode=mode, hx = 10.0 .* rand(8) .- 5, hy = 10.0 .* rand(8) .- 5,  dip=dip, azimuth=azimuth, taperz=(0.3,0.3))
+    lhs, rhs = dot_product_test(A,rand(domain(A)),rand(range(A)))
+    @test isapprox(lhs,rhs,rtol=1e-4)
+end
+
+@testset "JetSlantStackShiftSum3D, time-depth parity" begin
+    theta = collect(-45:1.0:45)
+    phi = [0.0, 90.0]
+    
+    # the ray parameters that would give the same results as theta/phi are
+    px = @. -tan(deg2rad(theta))
+    py = [0.0]
+    
+    Ad = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="depth", theta=theta, phi = phi, px=px, py=py, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64))
+    At1 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="time" , theta=theta, phi = phi, px=px, py=py, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64))
+    At2 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="time" , theta=theta, phi = phi, px=py, py=px, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64))
+    m = zeros(domain(Ad))
+    for i = 1:129
+        m[div(i,4)+20,i] = 1.0
+    end
+    m[16,64] = 1
+    m[32,96] = 1
+    dd = Ad*m
+    dt1 = At1*m
+    dt2 = At2*m
+
+    err_d1 = maximum(abs.(dd[:,:,1] .- dt1[:,:,1])) / maximum(abs.(dd[:,:,1]))
+    err_d2 = maximum(abs.(dd[:,:,2] .- dt2[:,1,:])) / maximum(abs.(dd[:,:,2]))
+
+    @show err_d1, err_d2
+    @test err_d1 < 1e-7
+    @test err_d2 < 1e-7
+end
+
+@testset "JetSlantStackShiftSum3D, dip parity" begin
+    theta = collect(-45:1.0:45)
+    phi = collect(0.0:45.0:135.0)
+    
+    A1 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="depth", theta=theta, phi = phi, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64), dip = 0.0, azimuth = 0.0)
+    A2 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="depth", theta=theta, phi = phi, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64), dip = zeros(64), azimuth = zeros(64))
+    m = zeros(domain(A1))
+    for i = 1:129
+        m[div(i,4)+20,i] = 1.0
+    end
+    m[16,64] = 1
+    m[32,96] = 1
+    d1 = A1*m
+    d2 = A2*m
+
+    err_d = maximum(abs.(d1 .- d2)) / maximum(abs.(d1))
+
+    @show err_d
+    @test err_d < 1e-7
+
+    A1 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="depth", theta=theta, phi = phi, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64), dip = 30.0, azimuth = 45.0)
+    A2 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, 129); mode="depth", theta=theta, phi = phi, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(-0.64:0.01:0.64), dip = 30 .+ zeros(64), azimuth = 45 .+ zeros(64))
+
+    d1 = A1*m
+    d2 = A2*m
+
+    err_d = maximum(abs.(d1 .- d2)) / maximum(abs.(d1))
+
+    @show err_d
+    @test err_d < 1e-7
+end
+
+@testset "JetSlantStackShiftSum3D, offsets parity" begin
+    theta = collect(-45:5.0:45)
+    phi = collect(0.0:45.0:135.0)
+
+    hx_reg = collect(-0.6:0.1:0.6)
+    hy_reg = collect(-0.3:0.1:0.2)
+    nhx = length(hx_reg)
+    nhy = length(hy_reg)
+    nh = nhx * nhy
+    hx_irreg = repeat(hx_reg, outer=nhy)
+    hy_irreg = repeat(hy_reg, inner=nhx)
+    
+    A1 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, nhx, nhy); mode="depth", theta=theta, phi = phi, dz=0.01, hx=hx_reg, hy=hy_reg, dip = 30.0, azimuth = 45.0)
+    A2 = JopSlantStackShiftSum3D(JetSpace(Float64, 64, nh); mode="depth", theta=theta, phi = phi, dz=0.01, hx=hx_irreg, hy=hy_irreg, dip = 30.0, azimuth = 45.0)
+    m1 = zeros(domain(A1))
+    for i = 1:nhy
+        m1[div(i,4)+20,:,i] .= 1.0
+    end
+    m1[16,3,2] = 1
+    m1[32,2, 4] = -1
+    m2 = reshape(m1, 64, nh)
+    
+    d1 = A1*m1
+    d2 = A2*m2
+
+    err_d = maximum(abs.(d1 .- d2)) / maximum(abs.(d1))
+
+    @show err_d
+    @test err_d < 1e-7
+end
+
+@testset "JetSlantStack3D vs JetSlantStackShiftSum3D, parity" for dip in (0.0, 10.0)
+    theta = collect(-45:1.0:45)
+    phi = collect(0.0:45.0:135.0)
+    azimuth = 45.0
+    A1 = JopSlantStack3D(JetSpace(Float64, 128, 129, 5); mode="depth", theta=theta, phi=phi, dz=0.01, hx0=-0.64, dhx=0.01, hy0=0, dhy=0.01, dip=dip, azimuth=azimuth)
+    A2 = JopSlantStackShiftSum3D(JetSpace(Float64, 128, 129, 5); mode="depth" , theta=theta, phi=phi, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(0:0.01:0.04), dip=dip, azimuth=azimuth)
+    m = zeros(domain(A1))
+    for i = 1:129
+        m[div(i,4)+20,i,1] = 1.0
+        m[div(i,3)+10,i,2] = -0.5
+    end
+    m[16,64,1] = 1
+    m[64,96,5] = -0.2
+    d1 = A1*m
+    d2 = A2*m
+
+    # check cosine similarity of the results since the two operators are not exactly the same
+    similarity = dot(d1, d2) / (norm(d1) * norm(d2))
+
+    @show similarity
+    @test similarity > 0.9
+end
+
+@testset "JopSlantStackShiftSum vs JopSlantStackShiftSum3D, parity" begin
+    theta = collect(-45:1.0:45)
+    phi = [0.0]
+    A1 = JopSlantStackShiftSum(JetSpace(Float64, 128, 129); mode="depth" , theta=theta, dz=0.01, h=collect(-0.64:0.01:0.64))
+    A2 = JopSlantStackShiftSum3D(JetSpace(Float64, 128, 129, 5); mode="depth" , theta=theta, phi=phi, dz=0.01, hx=collect(-0.64:0.01:0.64), hy=collect(0:0.01:0.04))
+    
+    m = zeros(domain(A2))
+    for i = 1:129
+        m[div(i,4)+20,i,1] = 1.0
+        m[div(i,3)+10,i,2] = -0.5
+    end
+    m[16,64,1] = 1
+    m[16,64,5] = -0.2
+    
+    d2 = A2*m
+    d1 = zeros(eltype(d2), size(d2)...)
+    for i = 1:size(m,3)
+        d1[:,:,1] .+= A1*view(m,:, :, i)
+    end
+
+    err_d = maximum(abs.(d1 .- d2)) / maximum(abs.(d1))
+    @show err_d
+    @test err_d < 1e-7
+end
\ No newline at end of file