From 9d8d246e9257add85b517b1516be28d58199789c Mon Sep 17 00:00:00 2001
From: Akshat Sharma <akshatsharma17112005@gmail.com>
Date: Sun, 1 Mar 2026 20:19:48 +0530
Subject: [PATCH] feat(transforms): add cosmological redshift and Doppler shift
 transforms; fix Unitful velocity conversion; add docs and tests (Issue #28)

---
 docs/src/transforms.md       |  40 ++++++++++
 src/Spectra.jl               |   1 +
 src/transforms/redshift.jl   | 105 +++++++++++++++++++++++++
 src/transforms/transforms.jl |   1 +
 test/transforms/redshift.jl  | 148 +++++++++++++++++++++++++++++++++++
 5 files changed, 295 insertions(+)
 create mode 100644 src/transforms/redshift.jl
 create mode 100644 test/transforms/redshift.jl

diff --git a/docs/src/transforms.md b/docs/src/transforms.md
index 571792a..091a718 100644
--- a/docs/src/transforms.md
+++ b/docs/src/transforms.md
@@ -60,6 +60,46 @@ deredden
 deredden!
 ```
 
+## Redshift
+
+Spectra can be shifted in wavelength via cosmological redshift or Doppler velocity. Both transformations operate only on the spectral axis; flux values are preserved as-is.
+
+**Cosmological redshift** shifts by a dimensionless parameter `z`:
+
+```jldoctest
+julia> using Spectra
+
+julia> spec = spectrum(collect(4000.0:1000.0:8000.0), ones(5));
+
+julia> shifted = redshift(spec, 0.5);
+
+julia> spectral_axis(shifted) ≈ spectral_axis(spec) .* 1.5
+true
+```
+
+**Doppler shift** shifts by a radial velocity `v`. Pass a `Unitful` velocity or a plain number (interpreted as m/s). Set `relativistic=true` for the full relativistic formula:
+
+```jldoctest
+julia> using Spectra, Unitful
+
+julia> spec = spectrum(collect(4000.0:1000.0:8000.0), ones(5));
+
+julia> shifted = doppler_shift(spec, 100u"km/s");
+
+julia> shifted_rel = doppler_shift(spec, 100u"km/s"; relativistic=true);
+```
+
+Both `redshift` and `doppler_shift` return a new spectrum. In-place variants `redshift!` and `doppler_shift!` are also available.
+
+### API/Reference
+
+```@docs
+redshift
+redshift!
+doppler_shift
+doppler_shift!
+```
+
 ## Resampling
 
 External interpolators, e.g., from [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) or [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl), can be used to resample spectra onto a given wavelength grid. Starting with a sample spectrum `spec`, we first create a [`SpectrumResampler`](@ref) object `resampler` which stores the initial spectrum and interpolator `interp` together. We then apply this object to the wavelength grid of our choice to produce the resampled spectrum. We show example usage in the docstring below:
diff --git a/src/Spectra.jl b/src/Spectra.jl
index 4c61508..d350584 100644
--- a/src/Spectra.jl
+++ b/src/Spectra.jl
@@ -8,6 +8,7 @@ export SingleSpectrum, IFUSpectrum, EchelleSpectrum
 
 # Transforms
 export SpectrumResampler, redden, redden!, deredden, deredden!
+export redshift, redshift!, doppler_shift, doppler_shift!
 
 # Utilities
 export blackbody #, line_flux, equivalent_width
diff --git a/src/transforms/redshift.jl b/src/transforms/redshift.jl
new file mode 100644
index 0000000..a499a8f
--- /dev/null
+++ b/src/transforms/redshift.jl
@@ -0,0 +1,105 @@
+const _c_mps = Float64(ustrip(c_0))
+
+_velocity_ratio(v::Real) = v / _c_mps
+_velocity_ratio(v::Quantity) = ustrip(u"m/s", v) / _c_mps
+
+"""
+    redshift!(spec::AbstractSpectrum, z::Real)
+
+In-place version of [`redshift`](@ref).
+"""
+function redshift!(spec::AbstractSpectrum, z::Real)
+    spec.spectral_axis = spectral_axis(spec) .* (1 + z)
+    return spec
+end
+
+"""
+    redshift(spec::AbstractSpectrum, z::Real)
+
+Apply a cosmological redshift to a spectrum, returning a new spectrum with shifted wavelengths.
+
+The observed wavelength is related to the rest-frame wavelength by
+
+``λ_\\mathrm{obs} = λ_\\mathrm{rest} \\cdot (1 + z)``
+
+where `z` is the cosmological redshift parameter. Only the spectral axis is transformed;
+flux density values are not corrected for the stretching of the wavelength bins.
+
+# Arguments
+- `spec`: The input spectrum.
+- `z`: Redshift parameter. Positive values redshift (longer wavelengths),
+  negative values blueshift (shorter wavelengths).
+
+# Examples
+```jldoctest
+julia> spec = spectrum(collect(4000.0:1000.0:8000.0), ones(5));
+
+julia> shifted = redshift(spec, 0.1);
+
+julia> spectral_axis(shifted) ≈ spectral_axis(spec) .* 1.1
+true
+```
+
+See also [`doppler_shift`](@ref) for velocity-based Doppler shifting.
+"""
+function redshift(spec::AbstractSpectrum, z::Real)
+    shifted = deepcopy(spec)
+    redshift!(shifted, z)
+    return shifted
+end
+
+"""
+    doppler_shift!(spec::AbstractSpectrum, v; relativistic=false)
+
+In-place version of [`doppler_shift`](@ref).
+"""
+function doppler_shift!(spec::AbstractSpectrum, v; relativistic::Bool = false)
+    β = _velocity_ratio(v)
+    if relativistic
+        spec.spectral_axis = spectral_axis(spec) .* sqrt((1 + β) / (1 - β))
+    else
+        spec.spectral_axis = spectral_axis(spec) .* (1 + β)
+    end
+    return spec
+end
+
+"""
+    doppler_shift(spec::AbstractSpectrum, v; relativistic=false)
+
+Apply a Doppler shift to a spectrum, returning a new spectrum with shifted wavelengths.
+
+**Non-relativistic** (default):
+
+``λ_\\mathrm{obs} = λ_\\mathrm{rest} \\cdot \\left(1 + \\frac{v}{c}\\right)``
+
+**Relativistic** (`relativistic=true`):
+
+``λ_\\mathrm{obs} = λ_\\mathrm{rest} \\cdot \\sqrt{\\frac{1 + v/c}{1 - v/c}}``
+
+Only the spectral axis is transformed; flux density values are preserved as-is.
+
+# Arguments
+- `spec`: The input spectrum.
+- `v`: Radial velocity. A `Unitful.Quantity` with velocity dimensions is automatically
+  converted using the speed of light in vacuum. A plain `Real` is interpreted as m/s.
+- `relativistic`: If `true`, use the relativistic Doppler formula. Default is `false`.
+
+# Examples
+```jldoctest
+julia> spec = spectrum(collect(4000.0:1000.0:8000.0), ones(5));
+
+julia> v = 1000.0;  # 1000 m/s
+
+julia> shifted = doppler_shift(spec, v);
+
+julia> spectral_axis(shifted) ≈ spectral_axis(spec) .* (1 + v / 299792458.0)
+true
+```
+
+See also [`redshift`](@ref) for cosmological redshift.
+"""
+function doppler_shift(spec::AbstractSpectrum, v; relativistic::Bool = false)
+    shifted = deepcopy(spec)
+    doppler_shift!(shifted, v; relativistic)
+    return shifted
+end
diff --git a/src/transforms/transforms.jl b/src/transforms/transforms.jl
index 91ed3e5..fcd974a 100644
--- a/src/transforms/transforms.jl
+++ b/src/transforms/transforms.jl
@@ -1,4 +1,5 @@
 include("redden.jl")
+include("redshift.jl")
 
 # ## Broadening ops
 # include("kernels.jl")
diff --git a/test/transforms/redshift.jl b/test/transforms/redshift.jl
new file mode 100644
index 0000000..42cea12
--- /dev/null
+++ b/test/transforms/redshift.jl
@@ -0,0 +1,148 @@
+using Spectra:
+    redshift,
+    redshift!,
+    doppler_shift,
+    doppler_shift!
+
+const C_MPS = 299792458.0
+
+function mock_spectrum_redshift(; use_units::Bool = false)
+    wave = collect(range(4000.0, 8000.0, length = 100))
+    flux = @. 1e14 / (wave^5 * (exp(1 / (wave * 6700)) - 1))
+    if use_units
+        wave *= u"angstrom"
+        flux *= u"erg/s/cm^2/angstrom"
+    end
+    spectrum(wave, flux, name = "Test Redshift Spectrum")
+end
+
+@testset "redshift" begin
+    @testset "Cosmological redshift z=$z" for z in [0.0, 0.1, 0.5, 1.0, 2.0]
+        spec = mock_spectrum_redshift()
+        original_wave = copy(spectral_axis(spec))
+        original_flux = copy(flux_axis(spec))
+        original_meta = deepcopy(spec.meta)
+
+        shifted = @inferred redshift(spec, z)
+
+        @test spectral_axis(shifted) ≈ original_wave .* (1 + z)
+        @test flux_axis(shifted) ≈ original_flux
+        @test shifted.meta == original_meta
+
+        # Original unchanged
+        @test spectral_axis(spec) ≈ original_wave
+        @test flux_axis(spec) ≈ original_flux
+    end
+
+    @testset "Cosmological redshift in-place" begin
+        spec = mock_spectrum_redshift()
+        expected_wave = spectral_axis(spec) .* 1.5
+
+        result = @inferred redshift!(spec, 0.5)
+        @test result === spec
+        @test spectral_axis(spec) ≈ expected_wave
+    end
+
+    @testset "Cosmological redshift roundtrip" begin
+        spec = mock_spectrum_redshift()
+        original_wave = copy(spectral_axis(spec))
+
+        shifted = redshift(spec, 1.0)
+        unshifted = redshift(shifted, -0.5)
+        @test spectral_axis(unshifted) ≈ original_wave
+    end
+
+    @testset "Cosmological redshift with Unitful" begin
+        spec = mock_spectrum_redshift(use_units = true)
+        original_wave = copy(spectral_axis(spec))
+
+        shifted = @inferred redshift(spec, 0.5)
+        @test spectral_axis(shifted) ≈ original_wave .* 1.5
+        @test flux_axis(shifted) ≈ flux_axis(spec)
+    end
+end
+
+@testset "doppler_shift" begin
+    @testset "Non-relativistic Doppler" begin
+        spec = mock_spectrum_redshift()
+        original_wave = copy(spectral_axis(spec))
+        original_flux = copy(flux_axis(spec))
+        original_meta = deepcopy(spec.meta)
+
+        v = 1e5  # 100 km/s in m/s
+        shifted = @inferred doppler_shift(spec, v)
+        expected = original_wave .* (1 + v / C_MPS)
+
+        @test spectral_axis(shifted) ≈ expected
+        @test flux_axis(shifted) ≈ original_flux
+        @test shifted.meta == original_meta
+
+        # Original unchanged
+        @test spectral_axis(spec) ≈ original_wave
+        @test flux_axis(spec) ≈ original_flux
+    end
+
+    @testset "Relativistic Doppler" begin
+        spec = mock_spectrum_redshift()
+        original_wave = copy(spectral_axis(spec))
+
+        v = 1e5  # 100 km/s in m/s
+        β = v / C_MPS
+        shifted = @inferred doppler_shift(spec, v; relativistic = true)
+        expected = original_wave .* sqrt((1 + β) / (1 - β))
+
+        @test spectral_axis(shifted) ≈ expected
+    end
+
+    @testset "Non-relativistic ≈ relativistic at low v" begin
+        spec = mock_spectrum_redshift()
+
+        v = 100.0  # 100 m/s — very non-relativistic
+        nr = doppler_shift(spec, v; relativistic = false)
+        rel = doppler_shift(spec, v; relativistic = true)
+        @test spectral_axis(nr) ≈ spectral_axis(rel) rtol = 1e-10
+    end
+
+    @testset "Doppler in-place" begin
+        spec = mock_spectrum_redshift()
+        v = 5e4
+        expected_wave = spectral_axis(spec) .* (1 + v / C_MPS)
+
+        result = @inferred doppler_shift!(spec, v)
+        @test result === spec
+        @test spectral_axis(spec) ≈ expected_wave
+    end
+
+    @testset "Doppler with Unitful velocity" begin
+        spec = mock_spectrum_redshift()
+        original_wave = copy(spectral_axis(spec))
+
+        v = 100.0u"km/s"
+        shifted = @inferred doppler_shift(spec, v)
+        β = 1e5 / C_MPS  # 100 km/s = 1e5 m/s
+        expected = original_wave .* (1 + β)
+        @test spectral_axis(shifted) ≈ expected
+
+        shifted_rel = @inferred doppler_shift(spec, v; relativistic = true)
+        expected_rel = original_wave .* sqrt((1 + β) / (1 - β))
+        @test spectral_axis(shifted_rel) ≈ expected_rel
+    end
+
+    @testset "Doppler with Unitful spectrum and velocity" begin
+        spec = mock_spectrum_redshift(use_units = true)
+        original_wave = copy(spectral_axis(spec))
+
+        v = 100.0u"km/s"
+        shifted = @inferred doppler_shift(spec, v)
+        β = 1e5 / C_MPS
+        @test spectral_axis(shifted) ≈ original_wave .* (1 + β)
+        @test flux_axis(shifted) ≈ flux_axis(spec)
+    end
+
+    @testset "Zero velocity" begin
+        spec = mock_spectrum_redshift()
+        shifted = doppler_shift(spec, 0.0)
+        @test spectral_axis(shifted) ≈ spectral_axis(spec)
+        @test flux_axis(shifted) ≈ flux_axis(spec)
+    end
+end