diff --git a/Project.toml b/Project.toml
index dc1a76f..681fc3a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,17 +1,31 @@
 name = "Gabs"
 uuid = "0eb812ee-a11f-4f5e-b8d4-bb8a44f06f50"
+<<<<<<< fix-issue-#59
+version = "1.3.4"
+authors = ["Andrew Kille"]
+=======
 authors = ["Andrew Kille"]
 version = "1.3.6"
+>>>>>>> main
 
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuantumInterface = "5717a53b-5d69-4fa3-b976-0bf2f97ca1e5"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SymplecticMatrices = "d07eab47-f98b-4620-8fe1-ee63e153d4ef"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
+TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 
 [weakdeps]
+<<<<<<< fix-issue-#59
+=======
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678"
+>>>>>>> main
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extensions]
@@ -20,14 +34,19 @@ QuantumOpticsBaseExt = "QuantumOpticsBase"
 StaticArraysExt = "StaticArrays"
 
 [compat]
+CairoMakie = "0.15.9"
 LinearAlgebra = "1.9"
 Makie = "0.21, 0.22, 0.23, 0.24"
+Plots = "1.41.6"
 QuantumInterface = "0.3.7, 0.3.8, 0.4.1"
 QuantumOpticsBase = "0.5"
 Random = "1.9"
 StaticArrays = "1.9.7"
 Symbolics = "6.27.0"
 SymplecticMatrices = "0.1.0"
+Test = "1.11.0"
+TestItemRunner = "1.1.4"
+TestItems = "1.0.0"
 julia = "1.6.7, 1.10.0"
 
 [extras]
diff --git a/src/channels.jl b/src/channels.jl
index 686fbba..71e7c95 100644
--- a/src/channels.jl
+++ b/src/channels.jl
@@ -220,6 +220,187 @@ function _amplifier(basis::QuadBlockBasis{N}, r::R, n::M) where {N<:Int,R<:Vecto
     return disp, transform, noise
 end
 
+""" 
+    classical_noise([Td=Vector{Float64}, Tt=Matrix{Float64},] basis::SymplecticBasis, n::Union{Int, Vector{Int}})
+    Gaussian channel describing the addition of classical noise to a Gaussian state. The channel is paramatrized
+    by noise parameter `n`, which can be a scalar (same noise on all modes) or a vector (different noise per mode).
+    The action is represented by the zero displacement vector and a noise matrix equal to `n*I'.
+
+    ## Example
+
+basis = QuadPairBasis(2)
+n_scalar = 0.8
+channel_scalar = classical_noise(basis, n_scalar)
+
+    SCALAR NOISE (n = 0.8)
+Displacement:
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+
+Transform Matrix:
+4×4 Matrix{Float64}:
+ 1.0  0.0  0.0  0.0
+ 0.0  1.0  0.0  0.0
+ 0.0  0.0  1.0  0.0
+ 0.0  0.0  0.0  1.0
+
+Noise Matrix:
+4×4 Matrix{Float64}:
+ 0.8  0.0  0.0  0.0
+ 0.0  0.8  0.0  0.0
+ 0.0  0.0  0.8  0.0
+ 0.0  0.0  0.0  0.8
+"""
+
+function classical_noise(::Type{Td}, ::Type{Tt}, basis::SymplecticBasis{N}, n::M; ħ = 2) where {Td,Tt,N<:Int,M}
+    disp, transform, noise = _classical_noise(basis, n)
+    return GaussianChannel(basis, Td(disp), Tt(transform), Tt(noise); ħ = ħ)
+end
+
+classical_noise(::Type{T}, basis::SymplecticBasis{N}, n::M; ħ = 2) where {T,N<:Int,M} = classical_noise(T, T, basis, n; ħ = ħ)
+
+function classical_noise(basis::SymplecticBasis{N}, n::M; ħ = 2) where {N<:Int,M}
+    disp, transform, noise = _classical_noise(basis, n)
+    return GaussianChannel(basis, disp, transform, noise; ħ = ħ)
+end
+
+# Internal function for scalar noise (same noise on all modes)
+function _classical_noise(basis::Union{QuadPairBasis{N},QuadBlockBasis{N}}, n::M) where {N<:Int,M}
+    nmodes = basis.nmodes
+    Rt = typeof(float(n))
+    disp = zeros(Rt, 2*nmodes) 
+    transform = Matrix{Rt}(I, 2*nmodes, 2*nmodes) # X = I
+    noise = Matrix{Rt}(n * I, 2*nmodes, 2*nmodes) # Y = n * I
+    return disp, transform, noise
+end
+
+# Internal function for vector noise (different noise per mode) 
+function _classical_noise(basis::QuadPairBasis{N}, n::M) where {N<:Int,M<:Vector}
+    nmodes = basis.nmodes
+    Mt = eltype(M)
+    disp = zeros(Mt, 2*nmodes) 
+    transform = Matrix{Mt}(I, 2*nmodes, 2*nmodes)
+    noise = zeros(Mt, 2*nmodes, 2*nmodes)
+    
+    @inbounds for i in Base.OneTo(nmodes)
+        ni = n[i]
+        noise[2*i-1, 2*i-1] = ni
+        noise[2*i, 2*i]     = ni
+    end
+    return disp, transform, noise
+end
+
+# Internal function for vector noise (different noise per mode)
+function _classical_noise(basis::QuadBlockBasis{N}, n::M) where {N<:Int,M<:Vector}
+    nmodes = basis.nmodeS
+    Mt = eltype(M)
+    disp = zeros(Mt, 2*nmodes) 
+    transform = Matrix{Mt}(I, 2*nmodes, 2*nmodes)
+    noise = zeros(Mt, 2*nmodes, 2*nmodes)
+    
+    @inbounds for i in Base.OneTo(nmodes)
+        ni = n[i]
+        noise[i, i]               = ni
+        noise[i+nmodes, i+nmodes] = ni
+    end
+    return disp, transform, noise
+end
+
+"""
+thermal_noise([Td=Vector{Float64}, Tt=Matrix{Float64},] basis::SymplecticBasis, η, c; ħ = 2)
+    Gaussian channel describing the interaction of a Gaussian state with a thermal environment. The channel is parameterized
+    by transmissivity `η` and mean thermal photon number `c`, which can be scalars (same for all modes) or vectors (different per mode).
+    The action is represented by a zero displacement vector, a transformation matrix equal to `sqrt(η)*I`, and a noise matrix equal to `(1-η)*(2c+1)*(ħ/2)*I`.
+    
+    ## Example
+    thermal_channel = thermal_noise(QuadPairBasis(2), 0.8, 0.6)
+    Displacement:
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 0.0
+ 0.0
+Transform Matrix:
+4×4 Matrix{Float64}:
+ 0.894427  0.0       0.0       0.0
+ 0.0       0.894427  0.0       0.0
+ 0.0       0.0       0.894427  0.0
+ 0.0       0.0       0.0       0.894427
+Noise Matrix:
+4×4 Matrix{Float64}:
+ 0.12  0.0   0.0   0.0
+ 0.0   0.12  0.0   0.0
+ 0.0   0.0   0.12  0.0
+ 0.0   0.0   0.0   0.12
+"""
+
+function thermal_noise(::Type{Td}, ::Type{Tt}, basis::SymplecticBasis{N}, η::M, c::L; ħ = 2) where {Td,Tt,N<:Int,M,L}
+    disp, transform, noise = _thermal_noise(basis, η, c)
+    return GaussianChannel(basis, Td(disp), Tt(transform), Tt(noise); ħ = ħ)
+end
+
+thermal_noise(::Type{T}, basis::SymplecticBasis{N}, η::M, c::L; ħ = 2) where {T,N<:Int,M,L} = thermal_noise(T, T, basis, η, c; ħ = ħ)
+
+function thermal_noise(basis::SymplecticBasis{N}, η::M, c::L; ħ = 2) where {N<:Int,M,L}
+    disp, transform, noise = _thermal_noise(basis, η, c)
+    return GaussianChannel(basis, disp, transform, noise; ħ = ħ)
+end
+
+# Internal function for scalar parameters 
+function _thermal_noise(basis::Union{QuadPairBasis{N},QuadBlockBasis{N}}, η::R, c::L) where {N<:Int,R,L}
+    nmodes = basis.nmodes
+    Rt = promote_type(typeof(float(η)), typeof(float(c))) 
+    disp = zeros(Rt, 2*nmodes) 
+    transform = Matrix{Rt}(sqrt(η) * I, 2*nmodes, 2*nmodes)
+    noise = Matrix{Rt}((1 - η) * c * I, 2*nmodes, 2*nmodes)
+    return disp, transform, noise
+end
+
+# Internal function for vector parameters 
+function _thermal_noise(basis::QuadPairBasis{N}, η::R, c::R) where {N<:Int,R<:Vector}
+    nmodes = basis.nmodes
+    Rt = eltype(R)
+    disp = zeros(Rt, 2*nmodes) 
+    transform = zeros(Rt, 2*nmodes, 2*nmodes)
+    noise = zeros(Rt, 2*nmodes, 2*nmodes)
+    
+    @inbounds for i in Base.OneTo(nmodes)
+        sqrt_η = sqrt(η[i])
+        noise_val = (1 - η[i]) * c[i]
+        
+        transform[2*i-1, 2*i-1] = sqrt_η
+        transform[2*i, 2*i]     = sqrt_η
+
+        noise[2*i-1, 2*i-1] = noise_val
+        noise[2*i, 2*i]     = noise_val
+    end
+    return disp, transform, noise
+end
+
+# Internal function for vector parameters 
+function _thermal_noise(basis::QuadBlockBasis{N}, η::R, c::R) where {N<:Int,R<:Vector}
+    nmodes = basis.nmodes
+    Rt = eltype(R)
+    disp = zeros(Rt, 2*nmodes) 
+    transform = zeros(Rt, 2*nmodes, 2*nmodes)
+    noise = zeros(Rt, 2*nmodes, 2*nmodes)
+    
+    @inbounds for i in Base.OneTo(nmodes)
+        sqrt_η = sqrt(η[i])
+        noise_val = (1 - η[i]) * c[i]
+
+        transform[i, i]               = sqrt_η
+        transform[i+nmodes, i+nmodes] = sqrt_η
+
+        noise[i, i]               = noise_val
+        noise[i+nmodes, i+nmodes] = noise_val
+    end
+    return disp, transform, noise
+end
+
 ##
 # Predefined operations on Gaussian channels
 ##
@@ -494,4 +675,5 @@ function changebasis(::Type{B1}, op::GaussianChannel{B2,D,S}) where {B1<:QuadPai
     return GaussianChannel(B1(nmodes), disp, transform, noise)
 end
 changebasis(::Type{<:QuadBlockBasis}, op::GaussianChannel{<:QuadBlockBasis,D,S}) where {D,S} = op
-changebasis(::Type{<:QuadPairBasis}, op::GaussianChannel{<:QuadPairBasis,D,S}) where {D,S} = op
\ No newline at end of file
+changebasis(::Type{<:QuadPairBasis}, op::GaussianChannel{<:QuadPairBasis,D,S}) where {D,S} = op
+
diff --git a/test/test_channels.jl b/test/test_channels.jl
index 9c4fe9c..6bdfb18 100644
--- a/test/test_channels.jl
+++ b/test/test_channels.jl
@@ -1,6 +1,7 @@
 @testitem "Channels" begin
     using Gabs
     using StaticArrays
+    import Gabs: classical_noise, thermal_noise
 
     nmodes = rand(1:5)
     qpairbasis = QuadPairBasis(nmodes)
@@ -199,4 +200,49 @@
 
         @test_throws AssertionError embed(QuadPairBasis(3), [1, 2, 3, 4], displace(QuadPairBasis(1), α, zeros(2, 2)))
     end
+
+
+@testset "classical noise channel" begin
+        n = rand(1:10)
+        ns = rand(1:10, nmodes)
+        
+        op_pair = classical_noise(qpairbasis, n)
+        op_block = classical_noise(qblockbasis, n)
+        
+        @test op_pair isa GaussianChannel && op_block isa GaussianChannel
+        @test classical_noise(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, qpairbasis, n) isa GaussianChannel
+        @test classical_noise(Array, qpairbasis, n) isa GaussianChannel
+        
+        @test op_pair == changebasis(QuadPairBasis, op_block) && op_block == changebasis(QuadBlockBasis, op_pair)
+        @test op_pair == changebasis(QuadPairBasis, op_pair) && op_block == changebasis(QuadBlockBasis, op_block)
+        
+        @test classical_noise(qblockbasis, n) == changebasis(QuadBlockBasis, op_pair)
+        @test classical_noise(qblockbasis, ns) == changebasis(QuadBlockBasis, classical_noise(qpairbasis, ns))
+        
+        @test isgaussian(op_pair, atol = 1e-4)
+        @test op_pair.ħ == 2 && op_block.ħ == 2
+    end
+
+    @testset "thermal noise channel" begin
+        eta = rand(Float64)
+        etas = rand(Float64, nmodes)
+        c = rand(1.0:10.0)
+        cs = rand(1.0:10.0, nmodes)
+        
+        op_pair = thermal_noise(qpairbasis, eta, c)
+        op_block = thermal_noise(qblockbasis, eta, c)
+        
+        @test op_pair isa GaussianChannel && op_block isa GaussianChannel
+        @test thermal_noise(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, qpairbasis, eta, c) isa GaussianChannel
+        @test thermal_noise(Array, qpairbasis, eta, c) isa GaussianChannel
+        
+        @test op_pair == changebasis(QuadPairBasis, op_block) && op_block == changebasis(QuadBlockBasis, op_pair)
+        @test op_pair == changebasis(QuadPairBasis, op_pair) && op_block == changebasis(QuadBlockBasis, op_block)
+        
+        @test thermal_noise(qblockbasis, eta, c) == changebasis(QuadBlockBasis, op_pair)
+        @test thermal_noise(qblockbasis, etas, cs) == changebasis(QuadBlockBasis, thermal_noise(qpairbasis, etas, cs))
+        
+        @test isgaussian(op_pair, atol = 1e-4)
+        @test op_pair.ħ == 2 && op_block.ħ == 2
+    end
 end
\ No newline at end of file