diff --git a/ext/ABCDMatrixOpticsPlotExt.jl b/ext/ABCDMatrixOpticsPlotExt.jl
index cfde7db..8fa482d 100644
--- a/ext/ABCDMatrixOpticsPlotExt.jl
+++ b/ext/ABCDMatrixOpticsPlotExt.jl
@@ -2,6 +2,12 @@ module ABCDMatrixOpticsPlotExt
 
 
 using ABCDMatrixOptics, Colors, Interpolations, RecipesBase
+export beamplot3d
+
+# mutable struct WithGaussianBeam3d
+	# system::Vector{<:Elements}
+	# beam::GaussianBeam3d
+# end
 
 mutable struct WithGaussianBeam
     system::Vector{<:Element}
@@ -14,6 +20,47 @@ mutable struct WithGeometricBeam
     beam::GeometricBeam
 end
 
+# did not succed to plot 2 shape with a single @receipe ..
+function beamplot3d(system::Vector{<:Elements}, beam::ABCDMatrixOptics.GaussianBeam3d)
+    ds = ABCDMatrixOptics2.discretize(system, 200)
+    N = length(ds) + 1
+    Tz = typeof(float(beam.zpos))
+    wsags = Vector{Tz}(undef, N)
+    wtans = Vector{Tz}(undef, N)
+    zs = Vector{Tz}(undef, N)
+    wsags[1] = ABCDMatrixOptics.wsag(beam)
+    wtans[1] = ABCDMatrixOptics.wtan(beam)
+    zs[1] = beam.zpos
+    beam = beam
+    for i = 1:length(ds)
+        beam = ds[i] * beam
+        wsags[i+1] = ABCDMatrixOptics.wsag(beam)
+        wtans[i+1] = ABCDMatrixOptics.wtan(beam)
+        zs[i+1] = beam.zpos 
+    end
+    xs, ys_sag, ys_tan = vcat(zs, reverse(zs)), vcat(wsags, (-1.0) .* reverse(wsags)), vcat(wtans, (-1.0) .* reverse(wtans))
+    shape_sag = [(xs[i], ys_sag[i], ) for i in 1:length(xs)]
+    shape_tan = [(xs[i], ys_tan[i], ) for i in 1:length(xs)]
+    index_iterfaces = findall(x -> !(x isa FreeSpace), ds)
+    # We look for all non-FreeSpace elements to display them in the final plot
+    z_interfaces = Vector{Tz}(undef, length(index_iterfaces))
+    for i=1:length(z_interfaces)
+        z_interfaces[i] = zs[index_iterfaces[i]]
+    end
+    plot(shape_sag, seriestype = :shape, 
+            linecolor = color(beam.λ),
+            fillcolor = color(beam.λ),
+            fillalpha = 0.1,
+            label = "Sag." * string(beam.λ),
+            xlabel = "Distance along Beam Axis",
+            yguide = "Extend along Beam Axis")
+    plot!(shape_tan, seriestype = :shape,
+            linecolor = color_tan(beam.λ),
+            fillcolor = color_tan(beam.λ),
+            fillalpha = 0.1,
+            label = "Tan." * string(beam.λ))
+    vline!(z_interfaces, primary = false, color = "black")
+end
 
 # type recipe, e.g. for `plot(WithBeam(system, beam))`
 @recipe f(::Type{WithGaussianBeam}, data::WithGaussianBeam) = begin
@@ -117,7 +164,23 @@ function color(λ)
     color = Colors.HSL(hue_of_λ(λ), 1.0, 0.5)
 end
 
+function color_tan(λ)
+	#slightly darker color for tangential display to be different
+    λmin = 380e-9
+    λmax = 800e-9
+    if λ < λmin
+        return Colors.HSL(260.0, 1.0, 0.15)
+    elseif λ > λmax
+        return Colors.HSL(0, 1.0, 0.15)
+    end
+    color = Colors.HSL(hue_of_λ(λ), 1.0, 0.40)
+end
+
 # user recipe, e.g. for `plot(system, beam)`
+# @recipe f(system::Vector{<:Elements}, beam::ABCDMatrixOptics.GaussianBeam3d) =
+	# WithGaussianBeam3d(system, beam)
+# @recipe f(system::Vector{<:Elements}, beam::ABCDMatrixOptics.GaussianBeam) =
+	# WithGaussianBeam3d(system, beam)
 @recipe f(system::Vector{<:Element}, beam::ABCDMatrixOptics.GaussianBeam) =
     WithGaussianBeam(system, beam)
 
diff --git a/src/beam.jl b/src/beam.jl
index c98dc7a..6850a87 100644
--- a/src/beam.jl
+++ b/src/beam.jl
@@ -1,4 +1,4 @@
-export AbstractBeam, GeometricBeam, GaussianBeam
+export AbstractBeam, GeometricBeam, GaussianBeam, GaussianBeam3d
 
 """
     abstract type AbstractBeam{T}
@@ -134,3 +134,113 @@ Return curvature `R` of the Gaussian beam.
 See also [`z`](@ref), [`zR`](@ref), [`w`](@ref), [`w0`](@ref).
 """
 R(beam::GaussianBeam) =  (z(beam) + zR(beam)^2 / z(beam)) 
+
+"""
+    struct GaussianBeam3d{T} <: AbstractBeam{T}
+        qsag
+		qtan
+        zpos
+        n
+        λ
+    end
+"""
+struct GaussianBeam3d{T} <:AbstractBeam{T}
+	qsag::Complex{T}
+	qtan::Complex{T}
+    zpos::T
+    n::T
+    λ::T
+end
+
+"""
+    GaussianBeam3d(w0=100e-6, w0sag = w0, w0tan = w0, z=0.0, zsag = z, ztan = z, n=1.0, λ=633e-9, zpos=0.0)
+
+Defines a physical Gaussian beam with different sagital and tangent plane.
+All parameters can be defined with keywords.
+
+Beam initial width `w0` and propagation distance `z` can be set for each `w0sag`, `zsag` and `w0tan`, `ztan` plane
+Refractive index `n` and vacuum wavelength `λ`.
+`zpos` defines a global position which is not used to calculate the q parameter but
+instead is for tracking purposes
+
+See also [`GaussianBeam`](@ref).
+
+
+## Example
+```jldoctest
+julia> GaussianBeam3d(w0=100e-6, z=12.0, n=1.0, λ=633e-9, zpos=0.0)
+GaussianBeam3d{Float64}(12.0 + 0.049630215696521214im, 12.0 + 0.049630215696521214im, 0.0, 1.0, 6.33e-7)
+```
+"""
+function GaussianBeam3d(; w0=100e-3, w0sag = w0, w0tan = w0, z=0.0, zsag = z, ztan = z, n=1.0, λ=633e-9, zpos=0.0)
+    w0sag, w0tan, zsag, ztan, n, λ, zpos = promote(w0sag, w0tan, zsag, ztan, n, λ, zpos)
+    return GaussianBeam3d{typeof(w0sag)}(zsag + 1im * (π * n * w0sag^2) / (λ), ztan + 1im * (π * n * w0tan^2) / (λ), zpos, n, λ)
+end
+
+"""
+    GaussianBeam3d(q; qsag = q, qtan = q, zpos, n=1.0, λ=633e-9)
+
+Returns a `GaussianBeam3d` defined by complex beam parameter `q`.
+
+
+## Example
+```jldoctest
+julia> GaussianBeam3d(12 + 1im * 1.0 * π * 100e-6^2 / 633e-9)
+GaussianBeam3d{Float64}(12.0 + 0.049630215696521214im, 12.0 + 0.049630215696521214im, 0.0, 1.0, 6.33e-7) 
+```
+"""
+function GaussianBeam3d(q; qsag = q, qtan = q, zpos=0.0, n=1.0, λ=633e-9)
+    zpos, n, λ = promote(zpos, n, λ)
+    qsag, qtan = Complex{typeof(zpos)}(qsag), Complex{typeof(zpos)}(qtan)
+    return GaussianBeam3d{real(typeof(qsag))}(qsag, qtan, zpos, n, λ)
+end
+
+"""
+    zsag(beam::GaussianBeam3d)
+
+Return the distance to the sagital beam waist `zsag` of the Gaussian beam. 
+All fsag where f is a function also exists for ftan
+
+"""
+zsag(beam::GaussianBeam3d) = real(beam.qsag)
+ztan(beam::GaussianBeam3d) = real(beam.qtan)
+
+"""
+    zRsag(beam::GaussianBeam3d)
+
+Return Rayleigh length `zRsag` in the sagital plane of the Gaussian beam.
+
+See also [`z`](@ref), [`w`](@ref), [`R`](@ref), [`w0`](@ref).
+"""
+zRsag(beam::GaussianBeam3d) = imag(beam.qsag)
+zRtan(beam::GaussianBeam3d) = imag(beam.qtan)
+
+"""
+    w0sag(beam::GaussianBeam3d)
+
+Return `w0` in the sagital plane of the Gaussian beam.
+
+See also [`z`](@ref), [`zR`](@ref), [`w`](@ref), [`R`](@ref).
+"""
+w0sag(beam::GaussianBeam3d) = sqrt(zRsag(beam) * beam.λ / (π * beam.n))
+w0tan(beam::GaussianBeam3d) = sqrt(zRtan(beam) * beam.λ / (π * beam.n))
+
+"""
+    wsag(beam::GaussianBeam3d)
+
+Return current width `w` in the sagital plane of the Gaussian beam.
+
+See also [`z`](@ref), [`zR`](@ref), [`R`](@ref), [`w0`](@ref).
+"""
+wsag(beam::GaussianBeam3d) = w0sag(beam) * sqrt(1 + (zsag(beam) / zRsag(beam))^2) 
+wtan(beam::GaussianBeam3d) = w0tan(beam) * sqrt(1 + (ztan(beam) / zRtan(beam))^2) 
+
+"""
+    Rsag(beam::GaussianBeam3d)
+
+Return curvature `R` in the sagital plane of the Gaussian beam.
+
+See also [`z`](@ref), [`zR`](@ref), [`w`](@ref), [`w0`](@ref).
+"""
+Rsag(beam::GaussianBeam3d) =  (zsag(beam) + zRsag(beam)^2 / zsag(beam)) 
+Rtan(beam::GaussianBeam3d) =  (ztan(beam) + zRtan(beam)^2 / ztan(beam)) 
diff --git a/src/elements.jl b/src/elements.jl
index 8cb1bcf..c37687a 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -1,7 +1,11 @@
+export Elements
 export Element, FreeSpace, Interface, ThinLens, ThickLens, Mirror
+export Element3d, Interface3d, ThinLens3d, ThickLens3d, Mirror3d
 export transfer_matrix
 
-abstract type Element{T} end
+abstract type Elements{T} end
+abstract type Element{T} <: Elements{T} end
+abstract type Element3d{T} <: Elements{T} end
 
 
 """
@@ -40,6 +44,32 @@ function ThickLens(; R1, R2, t, n_lens=1.5, n1=1.0, n2=1.0)
     return ThickLens{typeof(R1)}(R1, R2, t, n_lens, n1, n2)
 end
 
+"""
+    ThickLens3d(θ, R1, R2, t; n_lens=1.5 n1=1.0, n2=1.0)
+
+Construct a `θ`-tilted thick lens with the keywords:
+* `R1` radius of curvature of first surface
+* `R2` radius of curvature of second surface
+* `t`: thickness of lens
+* `n_lens=1.5` refractive index of lens
+* `n1=1`: refractive index of the medium of the incoming side
+* `n2=1`: refractive index of the medium of the exiting side
+
+"""
+struct ThickLens3d{T<:Number} <: Element3d{T}
+	θ::T
+	R1::T
+    R2::T
+    t::T
+    n_lens::T
+    n1::T
+    n2::T
+end
+function ThickLens3d(θ; R1, R2, t, n_lens=1.5, n1=1.0, n2=1.0)
+    θ, R1, R2, t, n_lens, n1, n2 = promote(θ, R1, R2, t, n_lens, n1, n2) 
+    return ThickLens3d{typeof(R1)}(θ, R1, R2, t, n_lens, n1, n2)
+end
+
 struct Interface{T<:Number} <: Element{T}
     n1::T
     n2::T
@@ -63,6 +93,20 @@ and `n2` on the new medium.
 """
 Interface(; n1, n2, R=Inf) = Interface{promote_type(typeof(n1), typeof(n2), typeof(R))}(promote(n1, n2, R)...)
 
+struct Interface3d{T<:Number} <: Element3d{T}
+	θ::T
+    n1::T
+    n2::T
+    R::T
+end
+
+"""
+    Interface3d(θ, n1=1, n2=1.5, R=Inf)
+
+Creates a θ (rad.) tilted interface with radius `R` and with refractive index `n1` on the entering side
+and `n2` on the new medium. `R`=Inf leads to a flat interface.
+"""
+Interface3d(θ ; n1=1, n2=1.5, R=Inf) = Interface3d{promote_type(typeof(θ), typeof(n1), typeof(n2), typeof(R))}(promote(θ, n1, n2, R)...)
 
 
 @with_kw_noshow struct ThinLens{T<:Number} <: Element{T}
@@ -85,6 +129,29 @@ The lens refractive index is `n_lens` and the outer refractive index is `n`.
 """
 ThinLens(R1, R2, n_lens=1.5, n=1.0) = ThinLens(inv((n_lens - n) / n * (1/R1 - 1/R2)))
 
+"""
+	ThinLens3d(θ, f)
+
+Creates a `θ`-tilted thin lens with focal length `f`.
+"""
+struct ThinLens3d{T<:Number} <: Element3d{T}
+    θ::T
+	f::T
+end
+
+"""
+    ThinLens3d(θ; R1, R2, n_lens = 1.5, n = 1)
+
+Creates a `θ`-tilted thin lens defined by the first radius of curvature `R1`, the second `R2`.
+The lens refractive index is `n_lens` and the outer refractive index is `n`.
+R>0 for concave surface
+"""
+function ThinLens3d(θ; R1, R2, n_lens=1.5, n=1.0) 
+	f = inv((n_lens - n) / n * (-1/R1 + 1/R2))
+	return ThinLens3d(θ, f)
+end
+ThinLens3d(θ, f) = ThinLens3d(promote(θ, f)...)
+
 """
     Mirror(R=Inf)
 
@@ -95,6 +162,16 @@ Per default `Inf`, so a flat mirror.
     R::T=Inf
 end
 
+"""
+    Mirror3d(θ, R=Inf)
+
+θ-tilted mirror with radius of curvature `R`.
+Per default `Inf`, so a flat mirror.
+"""
+@with_kw_noshow struct Mirror3d{T} <: Element{T}
+    θ::T
+	R::T=Inf
+end
 
 # definitions of dz
 """
@@ -109,6 +186,10 @@ dz(e::Mirror{T}) where T = zero(T)
 dz(e::ThickLens) = e.t
 dz(e::Matrix) = Inf
 
+dz(e::Interface3d{T}) where T = zero(T)
+dz(e::ThinLens3d{T}) where T = zero(T)
+dz(e::Mirror3d{T}) where T = zero(T)
+dz(e::ThickLens3d) = e.t
 
 
 
@@ -126,7 +207,21 @@ transfer_matrix(e::ThickLens) = transfer_matrix([Interface(n1=e.n1, n2=e.n_lens,
                                  FreeSpace(e.t), 
                                  Interface(n1=e.n_lens, n2=e.n2, R=e.R2)])
 
-
+"""
+    transfer_matrix(element::Element3d)
+	
+Returns both the Ray Transfer matrices (ABCDsag, ABCDtag) associated with the given, optical element
+"""
+# https://doi.org/10.1364/AO.26.000427 R>0 for concave surface, with e.θ2 asin(e.n1 * sin(e.θ) / e.n2) (Snell-Descartes law)
+transfer_matrix(e::Interface3d) = ([1 0 ; (e.n2*cos(asin(e.n1 * sin(e.θ) / e.n2)) - e.n1*cos(e.θ) )/(e.n2*e.R) e.n1/e.n2 ],
+				[cos(asin(e.n1 * sin(e.θ) / e.n2))/cos(e.θ) 0 ; ((e.n2*cos(asin(e.n1 * sin(e.θ) / e.n2)) - e.n1*cos(e.θ)) / (e.R * e.n2))  (e.n1*cos(e.θ))/(e.n2*cos(asin(e.n1 * sin(e.θ) / e.n2)))] )
+transfer_matrix(e::ThinLens3d) = ([1 0 ; -1/(e.f * cos(e.θ)) 1],
+				[1 0 ; -cos(e.θ)/e.f 1])
+transfer_matrix(e::Mirror3d) = ([1 0 ; -2*cos(e.θ)/e.R 1] ,
+				[1 0 ; -2/(e.R * cos(e.θ)) 1] )
+transfer_matrix(e::ThickLens3d) = transfer_matrix([Interface3d(e.θ, n1=e.n1, n2 = e.n_lens, R=e.R1),
+							FreeSpace(e.t),
+							Interface3d(e.θ, n1=e.n_lens, n2 = e.n2, R=e.R2)])
 
 """
     transfer_matrix(elements)
@@ -138,7 +233,27 @@ iteration) of optical elements.
 function transfer_matrix(elements::Vector{<:Element})
     return mapfoldr(transfer_matrix, (a,b) -> b * a, elements)
 end
+function transfer_matrix(elements::Vector{<:Elements})
+    return mapfoldr(transfer_matrix, (a,b) -> transfer_2matrices_3d(a,b), elements)
+	return
+end
 
+function transfer_2matrices_3d(A::Tuple, B::Tuple)
+	# A : (ABCDsag, ABCDtan) , B : (ABCDsag, ABCDtan)
+	return (B[1] * A[1] , B[2]* A[2] )
+end
+function transfer_2matrices_3d(A::Matrix, B::Tuple)
+	# A : ABCD , B : (ABCDsag, ABCDtan)
+	return (B[1] * A , B[2] * A )
+end
+function transfer_2matrices_3d(A::Tuple, B::Matrix)
+	# A : (ABCDsag, ABCDtan) , B : ABCD 
+	return (B * A[1] , B * A[2])
+end
+function transfer_2matrices_3d(A::Matrix, B::Matrix)
+	#Both are single ABCD matrix
+	return (B * A, B * A)
+end
 
 """
     discretize(e)
@@ -147,8 +262,8 @@ Discretizes the elements for plots. Nothing is done expect for FreeSpace, which
 """
 
 discretize(e::FreeSpace, N::Int) = fill(FreeSpace(e.dz/N), N)
-discretize(e::Element, N::Int) = e
-discretize(els::Vector{<:Element}, N::Int) = vcat(discretize.(els, Ref(N))...)
+discretize(e::Elements, N::Int) = e
+discretize(els::Vector{<:Elements}, N::Int) = vcat(discretize.(els, Ref(N))...)
 
 
 """
diff --git a/src/propagate.jl b/src/propagate.jl
index d022742..1c90b70 100644
--- a/src/propagate.jl
+++ b/src/propagate.jl
@@ -1,9 +1,10 @@
 export propagate, trace
 
 """
-    propagate(e::Union{Element, Vector{<:Element}}, b)
+    propagate(e::Union{Elements, Vector{<:Elements}}, b)
 
-Propagate a beam `b` either by a single element `e::Element` or `Vector{<:Element}`.
+Propagate a beam `b` either by a single element `e::Element` or `e::Element3d` or `Vector{<:Elements}`.
+If a Gaussian beams propagates through a 3D-element, it will become a 3D-Gaussian beam.
 
 Return is the final beam.
 Also available as `e * b`.
@@ -35,15 +36,48 @@ end
 _n(::Type{T}, e::Element, b, q_new) where T = GaussianBeam{T}(q_new, b.zpos + dz(e), b.n, b.λ) 
 _n(::Type{T}, e::Interface, b, q_new) where T = GaussianBeam{T}(q_new, b.zpos + dz(e), e.n2, b.λ) 
 
+function propagate(e::Element3d, b::GaussianBeam{T}) where T
+	ABCDsag, ABCDtan = transfer_matrix(e)
+	Asag, Bsag, Csag, Dsag = ABCDsag[1,1], ABCDsag[1,2], ABCDsag[2,1], ABCDsag[2,2]
+    qsag_new = (Asag * b.q + Bsag) / (Csag * b.q + Dsag)
+    Atan, Btan, Ctan, Dtan = ABCDtan[1,1], ABCDtan[1,2], ABCDtan[2,1], ABCDtan[2,2]
+    qtan_new = (Atan * b.q + Btan) / (Ctan * b.q + Dtan)
+    return _n(T, e, b, qsag_new, qtan_new)
+end
+
+function propagate(e::Element3d, b::GaussianBeam3d{T}) where T
+	ABCDsag, ABCDtan = transfer_matrix(e)
+	Asag, Bsag, Csag, Dsag = ABCDsag[1,1], ABCDsag[1,2], ABCDsag[2,1], ABCDsag[2,2]
+    qsag_new = (Asag * b.qsag + Bsag) / (Csag * b.qsag + Dsag)
+    Atan, Btan, Ctan, Dtan = ABCDtan[1,1], ABCDtan[1,2], ABCDtan[2,1], ABCDtan[2,2]
+    qtan_new = (Atan * b.qtan + Btan) / (Ctan * b.qtan + Dtan)
+    return _n(T, e, b, qsag_new, qtan_new)
+end
+
+function propagate(e::Element, b::GaussianBeam3d{T}) where T
+	ABCD = transfer_matrix(e)
+	A, B, C, D = ABCD[1,1], ABCD[1,2], ABCD[2,1], ABCD[2,2]
+    qsag_new = (A * b.qsag + B) / (C * b.qsag + D)
+	qtan_new = (A * b.qtan + B) / (C * b.qtan + D)
+    return _n(T, e, b, qsag_new, qtan_new)
+end
+
+_n(::Type{T}, e::Element3d, b::GaussianBeam, qsag_new, qtan_new) where T = GaussianBeam3d{T}(qsag_new, qtan_new, b.zpos + dz(e), b.n, b.λ)
+_n(::Type{T}, e::Element3d, b::GaussianBeam3d, qsag_new, qtan_new) where T = GaussianBeam3d{T}(qsag_new, qtan_new, b.zpos + dz(e), b.n, b.λ)
+_n(::Type{T}, e::Element, b::GaussianBeam3d, qsag_new, qtan_new) where T = GaussianBeam3d{T}(qsag_new, qtan_new, b.zpos + dz(e), b.n, b.λ)
+_n(::Type{T}, e::Interface, b::GaussianBeam3d, q_new) where T = GaussianBeam3d{T}(qsag_new, qtan_new, b.zpos + dz(e), e.n2, b.λ)
 
 
-function propagate(es::Vector{<:Element}, b::AbstractBeam)
+function propagate(es::Vector{<:Elements}, b::AbstractBeam)
     return reduce((a,b) -> propagate(b, a), es, init=b)
 end
 
-Base.:*(e::Union{Matrix, Element, Vector{<:Element}}, b::AbstractBeam) = propagate(e, b)
+Base.:*(e::Union{Matrix, Elements, Vector{<:Elements}}, b::AbstractBeam) = propagate(e, b)
 Base.:*(a::Element, b::Element) = transfer_matrix(a) * transfer_matrix(b)
+Base.:*(a::Elements, b::Elements) = transfer_matrix([a, b])
 Base.:*(a::Matrix, b::Element) = a * transfer_matrix(b)
+Base.:*(a::Matrix, b::Element3d) = a * transfer_matrix(b)
+Base.:*(a::Matrix, b::Tuple{Matrix, Matrix}) = (a*b[1], a*b[2])
 Base.:*(a::Vector{<:Element}, b::Vector) = transfer_matrix(a) * b