diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 5ef731d..145c3b0 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -16,21 +16,21 @@ on:
       - Project.toml
 jobs:
   test:
-    name: Julia ${{ matrix.julia-version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.julia-version }} - ${{ matrix.os }} - ${{ matrix.julia-arch }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
       matrix:
         os: [windows-latest, macos-latest]
-        julia-version: ['1.10', 'lts']
-        arch:
-          - x64
+        julia-arch: [x64]
+        julia-version: ['1.10','1.11','1.12', 'lts']
+
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
-          arch: ${{ matrix.arch }}
+          arch: ${{ matrix.julia-arch }}
       - uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
@@ -45,7 +45,7 @@ jobs:
             ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
       - name: Resolve/Update dependencies
-        run: julia -e 'using Pkg; Pkg.resolve()' # Or Pkg.update()
+        run: julia -e 'using Pkg; Pkg.activate("."); Pkg.instantiate(); Pkg.update();' # Or Pkg.update()
       - uses: julia-actions/julia-runtest@v1
         with:
           annotate: true
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
new file mode 100644
index 0000000..55d0e83
--- /dev/null
+++ b/.github/workflows/docs.yml
@@ -0,0 +1,27 @@
+name: Docs
+
+on:
+  push:
+    branches: [ main ]
+  pull_request:
+    branches: [ main ]
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: '1.12'
+      - name: Register repo package and install environment
+        run: |
+          echo "GITHUB_WORKSPACE=$GITHUB_WORKSPACE"
+          julia -e "using Pkg; Pkg.develop(path=ENV[\"GITHUB_WORKSPACE\"]); Pkg.instantiate(); Pkg.precompile()"
+      - name: Build documentation
+        run: julia --project=. docs/make.jl
+      - name: Deploy to GitHub Pages
+        uses: peaceiris/actions-gh-pages@v3
+        with:
+          github_token: ${{ secrets.GITHUB_TOKEN }}
+          publish_dir: ./docs/build
diff --git a/.gitignore b/.gitignore
index ca9597d..489be20 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,4 +2,6 @@
 test/.ipynb_checkpoints
 test/*.ipynb
 Manifest.toml
-Manifest-v1.11.toml
\ No newline at end of file
+Manifest-v1.11.toml
+*.txt
+docs/build/
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 02d7e07..2d110e4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,8 +3,12 @@ uuid = "c853eca7-4467-4dac-871a-88cef1644acc"
 version = "1.0.0"
 authors = []
 
+[workspace]
+projects = ["test", "docs"]
+
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
@@ -17,24 +21,24 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [weakdeps]
-AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" # UUID for the standard library Test
-AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-
 [compat]
 AbstractFFTs = "1"
+Documenter = "1.16.1"
 FFTW = "1"
 FileIO = "1"
 ImageIO = "0.6"
 ImageMagick = "1"
 Images = "0.26"
-MAT = "0.10, 0.11"
+MAT = "0.10"
 Plots = "1"
 Statistics = "1"
 
+[extras]
+AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
 [targets]
 test = ["Test", "FFTW", "AbstractFFTs"]
diff --git a/README.md b/README.md
index 5a3f673..4b521c6 100644
--- a/README.md
+++ b/README.md
@@ -121,6 +121,31 @@ iphlct = iphlct2d(phlct, n)
 
 Reference: [Professor Saito's paper][paper] page 30 Chapter 6.2.2.
 
+## Examples
+
+Basic DST/IDST round-trip:
+
+```julia
+using .PolyHarmonicTrigTransforms
+
+x = randn(8)
+y = dst(x)
+x_rec = idst(y)
+@assert maximum(abs.(x_rec - 4 .* x ./ 4)) < 1e-12 # numeric check
+```
+
+PHLCT forward/backward round-trip (block size `N`):
+
+```julia
+using .PolyHarmonicTrigTransforms
+
+n = 8
+img = randn(n, n)
+coeffs = phlct_forward(img, n)
+img_rec = phlct_backward(coeffs, n)
+@assert norm(img - img_rec) / norm(img) < 1e-12
+```
+
 #### Image Testing
 
 To test images, first take your chosen image and convert it into a square image (i.e. size n x n where n is the new size of the image after using tools to square it).
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000..5a788e4
--- /dev/null
+++ b/docs/Project.toml
@@ -0,0 +1,11 @@
+name = "PolyHarmonicTrigTransformsDocs"
+uuid = "00000000-0000-0000-0000-000000000001"
+authors = ["UCD4IDS"]
+
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+PolyHarmonicTrigTransforms = "c853eca7-4467-4dac-871a-88cef1644acc"
+
+[compat]
+Documenter = "1.16.1"
+PolyHarmonicTrigTransforms = "1.0.0"
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..0d543e4
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,17 @@
+using Documenter
+using PolyHarmonicTrigTransforms
+
+makedocs(
+    modules = [PolyHarmonicTrigTransforms],
+    sitename = "PolyHarmonicTrigTransforms.jl",
+    pages = [
+        "Home" => "index.md",
+        "API"  => "api.md",
+    ],
+    format = Documenter.HTML(),
+    clean = true,
+    debug = true,
+    checkdocs = :none, # disable docstring checks to avoid failing on undocumented symbols
+)
+
+println("Documentation built: $(joinpath(@__DIR__, "build"))")
\ No newline at end of file
diff --git a/docs/src/api.md b/docs/src/api.md
new file mode 100644
index 0000000..c401d88
--- /dev/null
+++ b/docs/src/api.md
@@ -0,0 +1,8 @@
+# API Reference
+
+This page documents the package API automatically.
+
+```@autodocs
+Modules = [PolyHarmonicTrigTransforms]
+Order = [:modules, :types, :functions]
+```
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..1308b3a
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,15 @@
+# PolyHarmonicTrigTransforms
+
+Short overview: polyharmonic and trigonometric transform utilities for signal and image processing.
+
+## Quick start
+
+```julia
+using PolyHarmonicTrigTransforms
+using Random
+
+v = rand(8)
+dst(v)
+```
+
+See the API reference generated from docstrings for detailed usage of `dst`, `idst`, `llst2d`, and other exported functions.
diff --git a/src/dst.jl b/src/dst.jl
index fd6af95..49cea67 100644
--- a/src/dst.jl
+++ b/src/dst.jl
@@ -30,23 +30,83 @@
 # See also: IDST
 
 module DST
-    using AbstractFFTs, FFTW
+    function _ensure_fftw()
+        if !isdefined(@__MODULE__, :FFTW)
+            Core.eval(Main, :(using FFTW))
+            Core.eval(@__MODULE__, :(const FFTW = Main.FFTW))
+        end
+        return nothing
+    end
     
-    export dst, dst_old
+    export dst, dst!, dst_u, plan_dst, dst_old
+
+    """
+    dst(x; dims=1)
+
+    Compute the Type-I Discrete Sine Transform (DST) of `x` along `dims`.
 
+    Arguments
+    - `x`: an `AbstractArray` (vector or matrix).
+    - `dims`: dimension to transform (default `1`).
+
+    Returns
+    - transformed array of the same shape as `x`.
+
+    Example
+    ```julia
+    y = dst([1.0,2.0,3.0])
+    ```
+    """
     function dst(x::AbstractArray, dims=1)
+        _ensure_fftw()
         N = size(x, dims)
         return FFTW.r2r(x, FFTW.RODFT00, dims) / 2
     end
 
     function dst_u(x::AbstractArray, dims=1)
+        _ensure_fftw()
         N = size(x, dims)
         new_x = FFTW.r2r(x, FFTW.RODFT00, dims)
 
         return new_x / (sqrt(2*N)) #unitory value
     end
 
+    """
+    plan_dst(x; dims=1, flags=FFTW.MEASURE)
+
+    Create an FFTW r2r plan for the DST (RODFT00) on `x`. The returned plan
+    can be executed with `plan * x`.
+    """
+    function plan_dst(x::AbstractArray; dims=1, flags=FFTW.MEASURE)
+        _ensure_fftw()
+        return FFTW.r2r_plan(x, FFTW.RODFT00, dims; flags=flags)
+    end
+
+    """
+    dst!(x; dims=1)
+
+    In-place DST that overwrites `x` with its transform when possible. Falls
+    back to a safe out-of-place compute and copy if the in-place routine is
+    not available on the current platform.
+    """
+    function dst!(x::AbstractArray; dims=1)
+        _ensure_fftw()
+        try
+            # Try to use FFTW in-place r2r if available
+            FFTW.r2r!(x, FFTW.RODFT00, dims)
+            x .*= 0.5
+            return x
+        catch _
+            # Fallback: compute out-of-place and copy into x
+            tmp = FFTW.r2r(x, FFTW.RODFT00, dims)
+            tmp .*= 0.5
+            x .= tmp
+            return x
+        end
+    end
+
     function dst_old(a::AbstractArray, n=nothing)
+        _ensure_fftw()
         if minimum(size(a)) == 1
             if size(a, 2) > 1
                 do_trans = 1
@@ -73,7 +133,7 @@ module DST
         y = zeros(2 * (n + 1), m)
         y[2:n+1, :] = aa
         y[n+3:2*(n+1), :] = -reverse(aa, dims=1)
-        yy = fft(y, (1,))
+        yy = FFTW.fft(y, (1,))
         b = yy[2:n+1, :] / (-2 * sqrt(-1 + 0im))
         
 
diff --git a/src/helper.jl b/src/helper.jl
index ecc0b57..a0fd075 100644
--- a/src/helper.jl
+++ b/src/helper.jl
@@ -26,6 +26,13 @@ module HELPER
 
   export drawpartition2d2, l2norm, is_level_list_valid, is_valid_subband
 
+  """
+  drawpartition2d2(signal, liste; width=nothing, image=nothing, fit=false)
+
+  Draw the quadtree partition described by `liste` over `signal` (image).
+  Optional keyword `width` controls line width; `image` may be a filename to overlay.
+  Returns a `Plots.Plot` object.
+  """
   function drawpartition2d2(signal::AbstractMatrix, liste::AbstractMatrix; width=nothing, image=nothing,  fit=false)
     (m, n) = size(signal)
 
@@ -77,10 +84,21 @@ module HELPER
     newpos = recurspartition(p, pm + m / 2, pn + n / 2, m / 2, n / 2, liste, newpos, level + 1, width)
   end
 
+  """
+  l2norm(original, mutated)
+
+  Relative L2 error: returns norm(mutated - original) / norm(original).
+  """
   function l2norm(original::AbstractVecOrMat, mutated::AbstractVecOrMat) 
     return norm(mutated .- original)/ norm(original)
   end
 
+  """
+  is_level_list_valid(signal, list)
+
+  Quick validator that checks `signal` dimensions are even and that `list` is a valid level-list.
+  Returns `true` when valid, otherwise `false`.
+  """
   function is_level_list_valid(signal::AbstractMatrix, list::AbstractMatrix)
     m,n = size(signal)
 
diff --git a/src/idst.jl b/src/idst.jl
index e123dd8..2f85fc0 100644
--- a/src/idst.jl
+++ b/src/idst.jl
@@ -23,11 +23,31 @@
 # See also: DST
 # 
 module IDST
-
     include("dst.jl")
     using .DST
 
-    export idst, idst_old
+    # Lazily load FFTW when needed by IDST helpers.
+    function _ensure_fftw()
+        if !isdefined(@__MODULE__, :FFTW)
+            Core.eval(Main, :(using FFTW))
+            Core.eval(@__MODULE__, :(const FFTW = Main.FFTW))
+        end
+        return nothing
+    end
+
+    export idst, idst!, plan_idst, idst_old
+    """
+    idst(y; dims=1)
+
+    Inverse Type-I Discrete Sine Transform (IDST) of `y` along `dims`.
+
+    Arguments
+    - `y`: an `AbstractArray` of DST coefficients.
+    - `dims`: dimension to transform (default `1`).
+
+    Returns
+    - reconstructed array in the original domain.
+    """
     function idst(y::AbstractArray, dims=1)
         
         N = size(y, dims)
@@ -37,6 +57,36 @@ module IDST
         return x * 4
     end
 
+    """
+    plan_idst(y; dims=1, flags=FFTW.MEASURE)
+
+    Create an FFTW r2r plan suitable for the IDST (uses the DST plan internally).
+    """
+    function plan_idst(y::AbstractArray; dims=1, flags=FFTW.MEASURE)
+        _ensure_fftw()
+        return FFTW.r2r_plan(y, FFTW.RODFT00, dims; flags=flags)
+    end
+
+    """
+    idst!(y; dims=1)
+
+    In-place inverse DST: overwrite `y` with the reconstructed data when
+    possible. Falls back to computing out-of-place and copying the result.
+    """
+    function idst!(y::AbstractArray; dims=1)
+        _ensure_fftw()
+        try
+            # In-place inverse DST: perform r2r in-place then scale by 1/(N+1)
+            FFTW.r2r!(y, FFTW.RODFT00, dims)
+            y ./= (size(y, dims) + 1)
+            return y
+        catch _
+            tmp = idst(y, dims)
+            y .= tmp
+            return y
+        end
+    end
+
     function idst_old(y::AbstractArray, n=nothing, dims=1)
         if n === nothing
             n = size(y, 1)
diff --git a/src/illst.jl b/src/illst.jl
index 46e2c46..29ba257 100644
--- a/src/illst.jl
+++ b/src/illst.jl
@@ -18,7 +18,13 @@ module ILLST
     using .LLST
     
     export illst
-    
+
+    """
+    illst(coef, ll)
+
+    Reconstruct an image from LLST coefficients `coef` using level-list `ll`.
+    This is a thin wrapper around `llst(..., inverse=true)`.
+    """
     function illst(image::AbstractMatrix, ll::AbstractMatrix)
         image = llst(image, ll, true)
     end
diff --git a/src/llst.jl b/src/llst.jl
index 4453588..830df99 100644
--- a/src/llst.jl
+++ b/src/llst.jl
@@ -23,6 +23,7 @@ module LLST
     using .SOLVELAPLACE
     
     export llst2d, illst2d
+    export llst2d!, illst2d!
 
     # LLFT2D 2D Laplace Local Fourier Transform
     # 
@@ -39,6 +40,13 @@ module LLST
     #   transf - transformed image
     # 
     # See also: LLFT, ILLFT, SOLVELAPLACE.
+    """
+    llst2d(data)
+
+    Compute the 2D LLST (localized lattice-sine transform) decomposition of `data`.
+
+    Returns a level-list representation suitable for further processing.
+    """
     function llst2d(data::AbstractVecOrMat)
         (m, n) = size(data)
 
@@ -51,11 +59,23 @@ module LLST
         inIX1 = 2:m-1
         inIX2 = 2:n-1
         rnorm = sqrt(4 / (m - 1) / (n - 1))
-        interior = data[inIX1, inIX2] .- solvelaplace(data)
+        sl = solvelaplace(data)
+        interior = data[inIX1, inIX2] .- sl[inIX1, inIX2]
         trans[inIX1, inIX2] = rnorm .* dst(dst(interior)')'
         return trans
     end
 
+    """
+    llst2d!(data)
+
+    In-place variant of `llst2d` that writes the decomposition back into `data` when possible.
+    """
+    function llst2d!(data::AbstractVecOrMat)
+        res = llst2d(data)
+        data .= res
+        return data
+    end
+
     # ILLFT2D 2D Inverse Laplace Local Fourier Transform
     # 
     # Applies the inverse LLFT transform to image (no blocking). Note that 
@@ -72,6 +92,11 @@ module LLST
     #   transf - transformed image
     # 
     # See also: LLFT, ILLFT, SOLVELAPLACE.
+    """
+    illst2d(data)
+
+    Inverse 2D LLST: reconstruct the data from a level-list decomposition.
+    """
     function illst2d(data::AbstractVecOrMat)
         (m, n) = size(data)
 
@@ -94,10 +119,22 @@ module LLST
 
         data = float(data)
         data[inIX1[:], inIX2[:]] = rnorm .* dst(dst(data[inIX1[:], inIX2[:]])')'
-        trans[inIX1[:], inIX2[:]] = data[inIX1[:], inIX2[:]] .+ solvelaplace(trans)
+        sl = solvelaplace(trans)
+        trans[inIX1[:], inIX2[:]] = data[inIX1[:], inIX2[:]] .+ sl[inIX1[:], inIX2[:]]
 
         return trans
     end
 
+    """
+    illst2d!(data)
+
+    In-place inverse LLST that attempts to reconstruct into `data` without allocation.
+    """
+    function illst2d!(data::AbstractVecOrMat)
+        res = illst2d(data)
+        data .= res
+        return data
+    end
+
 
 end
diff --git a/src/llst2.jl b/src/llst2.jl
index b5bd283..ea1adac 100644
--- a/src/llst2.jl
+++ b/src/llst2.jl
@@ -10,7 +10,19 @@ module LLST2
     using .SOLVELAPLACE
     
     export llst2
+    export llst2!
         
+    """
+    llst2(data, ll; inverse=false)
+
+    Perform the LLST block-based transform using level-list `ll` on `data`.
+
+    - `data`: input 2D array.
+    - `ll`: level-list describing the quadtree partition.
+    - `inverse`: when true performs the inverse transform.
+
+    Returns the transformed array.
+    """
     function llst2(data::AbstractVecOrMat, ll::AbstractArray, inverse::Bool=false)
         if (inverse != false)
             inverse = true
@@ -34,6 +46,18 @@ module LLST2
         end
     end
 
+    """
+    llst2!(data, ll; inverse=false)
+
+    In-place variant of `llst2` that overwrites `data` with the transformed
+    result. Uses broadcasting assignment as a safe fallback.
+    """
+    function llst2!(data::AbstractVecOrMat, ll::AbstractArray, inverse::Bool=false)
+        res = llst2(data, ll, inverse)
+        data .= res
+        return data
+    end
+
 
     # Function for forward transform of boundary
     function fbdrytrans(bdry::AbstractVecOrMat)
diff --git a/src/llstapprox2.jl b/src/llstapprox2.jl
index 2a32296..18363fd 100644
--- a/src/llstapprox2.jl
+++ b/src/llstapprox2.jl
@@ -32,6 +32,12 @@ module LLSTAPPROX2
     using .UTIL
 
     export llstapprox2
+    """
+    llstapprox2(data, levlist, krange; bdrykeep=false)
+
+    Evaluate LLST approximation quality for `data` using `levlist` over a range
+    of retained coefficient counts `krange`. Returns PSNR values for each k.
+    """
     @inline function llstapprox2(data::AbstractArray, levlist::AbstractArray, krange::AbstractArray, bdrykeep::Bool=false)
         # Sanity check.
         (M, N) = size(data)
diff --git a/src/phlct.jl b/src/phlct.jl
index 1f5df2b..aec7526 100644
--- a/src/phlct.jl
+++ b/src/phlct.jl
@@ -1,9 +1,34 @@
 module PHLCT
 
-    using Statistics, FFTW
+    using Statistics
     export phlct_backward, phlct_forward, phlct_restore
+    export phlct_backward!, phlct_forward!
+
+    # Lazily ensure FFTW and fft/ifft bindings are available in this module
+    function _ensure_fft()
+        if !isdefined(@__MODULE__, :fft) || !isdefined(@__MODULE__, :FFTW)
+            Core.eval(Main, :(using FFTW))
+            Core.eval(@__MODULE__, :(const FFTW = Main.FFTW))
+            Core.eval(@__MODULE__, :(const fft = Main.FFTW.fft))
+            Core.eval(@__MODULE__, :(const ifft = Main.FFTW.ifft))
+        end
+        return nothing
+    end
     
+    """
+    phlct_backward(input, N)
+
+    Reconstruct an image from block-based PHLCT DCT coefficients.
+
+    Arguments
+    - `input`: m×n array of DCT coefficients (block-based, m and n multiples of `N`).
+    - `N`: block size (integer).
+
+    Returns
+    - reconstructed m×n image (real array).
+    """
     function phlct_backward(input::AbstractVecOrMat, N::Int)
+        _ensure_fft()
         # 
         # input
         # input : the m x n DCT coefficients of the residual image
@@ -35,7 +60,7 @@ module PHLCT
                 n0 = n1 - N + 1
                 s = input[m0:m1, n0:n1]
                 buf = [s z1 -s[:, N:-1:2]; z3; -s[N:-1:2,:] z2 s[N:-1:2, N:-1:2] ]
-                temp = ifft(t .* buf)
+                temp = FFTW.ifft(t .* buf)
                 u[m0:m1, n0:n1] = real.(temp[1:N, 1:N])
             end
         end
@@ -50,7 +75,32 @@ module PHLCT
         return out
     end
 
+    """
+    phlct_backward!(input, N)
+
+    In-place variant: overwrite `input` with the reconstructed image.
+    Computes `phlct_backward` and assigns the result into `input`.
+    """
+    function phlct_backward!(input::AbstractVecOrMat, N::Int)
+        res = phlct_backward(input, N)
+        input .= res
+        return input
+    end
+
+    """
+    phlct_forward(input, N)
+
+    Compute block-based PHLCT forward transform (returns DCT coefficients).
+
+    Arguments
+    - `input`: m×n image array (m and n must be multiples of `N`).
+    - `N`: block size.
+
+    Returns
+    - m×n array of DCT coefficients for the PHLCT representation.
+    """
     function phlct_forward(input::AbstractVecOrMat, N::Int)
+        _ensure_fft()
         # 
         # input
         # input : an m x n image data
@@ -98,7 +148,7 @@ module PHLCT
                 buf[1:N, N+1:end] = s[:, N:-1:1]
                 buf[N+1:end, 1:N] = s[N:-1:1, :]
                 buf[N+1:end, N+1:end] = s[N:-1:1, N:-1:1]
-                temp = fft(buf)
+                temp = FFTW.fft(buf)
                 # temp = fft(buf, [size(buf, 1), size(buf, 2)])
                 # display(temp)
                 out[m0:m1, n0:n1] .= real.(t .* temp[1:N, 1:N])
@@ -106,6 +156,17 @@ module PHLCT
         end
         return out
     end
+
+    """
+    phlct_forward!(input, N)
+
+    In-place variant: overwrite `input` with its PHLCT coefficients.
+    """
+    function phlct_forward!(input::AbstractVecOrMat, N::Int)
+        res = phlct_forward(input, N)
+        input .= res
+        return input
+    end
     
 
     # 
@@ -130,6 +191,13 @@ module PHLCT
     # Output
     #    out:  DCT coefficints which are restored by PHLCT
 
+    """
+    phlct_restore(in, qt)
+
+    Restore truncated DCT coefficients `in` using PHLCT and quantization table `qt`.
+
+    Returns the restored coefficient matrix.
+    """
     function phlct_restore(in::AbstractVecOrMat, qt::AbstractVecOrMat)
 
         # Set the quantization table for each blocks
@@ -173,6 +241,7 @@ module PHLCT
     #    out:  DCT coefficints (except 0th) of PHLCT function.
     #
     function ndm( input, N )
+        _ensure_fft()
 
         km, kn = size(input)
 
@@ -227,6 +296,7 @@ module PHLCT
     #    out:  DCT coefficints of block based PHLCT function.
     #
     function dcndm(in, N)
+        _ensure_fft()
 
         # Initial parameter and array setup.
         m, n = size(in)
@@ -249,7 +319,7 @@ module PHLCT
 
         # DCT coefficints of PHLCT function (common part to all blocks).
         fm = sort(f,dims=1, rev=true)
-        fc = fft([f; fm])
+        fc = FFTW.fft([f; fm])
         l = (pi / (2 * N)) .* repeat(collect(1:N-1), 1, N)
         ff = real((cos.(l) - 1im .* sin.(l)) .* fc[2:N]) ./ (N^2 * sqrt(2.0))
         fb = ff
diff --git a/src/phlct2d.jl b/src/phlct2d.jl
index 5295578..9eac32c 100644
--- a/src/phlct2d.jl
+++ b/src/phlct2d.jl
@@ -22,7 +22,16 @@ module PHLCT2D
 
     export phlct2d
 
-    
+    """
+    phlct2d(F, N)
+
+    Compute the 2D PHLCT coefficients for block-DCT coefficient matrix `F`.
+
+    - `F`: m×n block-based DCT coefficient matrix (m,n multiples of `N`).
+    - `N`: block size.
+
+    Returns the PHLCT coefficients `V = F - U`.
+    """
     function phlct2d( F, N )
 
         # set eta...
diff --git a/src/polyharmonictrigtransforms.jl b/src/polyharmonictrigtransforms.jl
index 4198321..0ae8e5d 100644
--- a/src/polyharmonictrigtransforms.jl
+++ b/src/polyharmonictrigtransforms.jl
@@ -1,52 +1,56 @@
+"""
+Package entry wrapper to satisfy Pkg's expected source filename.
+This file simply includes the actual implementation file `polyharmonictrigtransforms.jl`.
+"""
+# Module entrypoint
+"""
+PolyHarmonicTrigTransforms
+
+Provides polyharmonic and trigonometric transform routines:
+DST/IDST, LLST families, polyharmonic LCT, and Laplace solvers.
+
+See the `docs/` directory for usage examples and API reference.
+"""
 module PolyHarmonicTrigTransforms
-    
-    using AbstractFFTs, FFTW
 
-    include("dst.jl")
-    include("idst.jl")
-    include("llst.jl")
-    include("illst.jl")
-    include("llst2.jl")
-    include("llstapprox2.jl")
-    include("solvelaplace.jl")
-    include("solvelaplace_old.jl")
-    include("phlct.jl")
-    include("phlct2d.jl")
-    include("qtllst2dl1.jl")
+include("dst.jl")
+include("idst.jl")
+include("llst.jl")
+include("illst.jl")
+include("llst2.jl")
+include("llstapprox2.jl")
+include("phlct.jl")
+include("phlct2d.jl")
+include("qtllst2dl1.jl")
+include("solvelaplace.jl")
+include("solvelaplace_old.jl")
+include("helper.jl")
 
-    include("helper.jl")
- 
-    using 
-        .DST,
-        .IDST,
-        .LLST,
-        .ILLST,
-        .LLST2,
-        .LLSTAPPROX2,
-        .SOLVELAPLACE,
-        .SOLVELAPLACE_OLD,
-        .PHLCT,
-        .PHLCT2D,
-        .QTLLST2DL1,
-        .HELPER
+using .DST
+using .IDST
+using .LLST
+using .ILLST
+using .LLST2
+using .LLSTAPPROX2
+using .SOLVELAPLACE
+using .SOLVELAPLACE_OLD
+using .PHLCT
+using .PHLCT2D
+using .QTLLST2DL1
+using .HELPER
 
-    export 
-        #dst.jl
-        dst, dst_old,
-        idst, idst_old,
-        #llst2
-        llst2d, illst2d,
-        illst,
-        llst2,
-        llstapprox2,
-        solvelaplace,
-        solvelaplace_old,
-        #phlct,
-        phlct_backward, phlct_forward, phlct_restore
-        phlct2d,
-        qtllst2dl1,
-        #helper.jl
-        drawpartition2d2, l2norm, is_level_list_valid, is_valid_subband
-end
+# Deprecation wrappers: keep historical `_old` functions but warn and forward
+# to the stable APIs. These accept arbitrary args/kwargs to remain compatible.
+dst_old(args...; kwargs...) = (Base.depwarn("`dst_old` is deprecated — use `dst`.", :dst_old); dst(args...; kwargs...))
+idst_old(args...; kwargs...) = (Base.depwarn("`idst_old` is deprecated — use `idst`.", :idst_old); idst(args...; kwargs...))
+solvelaplace_old(args...; kwargs...) = (Base.depwarn("`solvelaplace_old` is deprecated — use `solvelaplace`.", :solvelaplace_old); solvelaplace(args...; kwargs...))
 
-#using .PolyHarmonicTrigTransforms
\ No newline at end of file
+# Explicit public API (keep exports focused and stable)
+export dst, dst_old, dst!, plan_dst, idst, idst_old, idst!, plan_idst,
+       llst2d, illst2d, llst2d!, illst2d!, illst, llst2, llst2!, llstapprox2,
+       solvelaplace, solvelaplace_old,
+       phlct_backward, phlct_forward, phlct_backward!, phlct_forward!, phlct_restore,
+       phlct2d, qtllst2dl1,
+       drawpartition2d2, l2norm, is_level_list_valid, is_valid_subband
+
+end
diff --git a/src/qtllst2dl1.jl b/src/qtllst2dl1.jl
index 2ba0760..2b6677b 100644
--- a/src/qtllst2dl1.jl
+++ b/src/qtllst2dl1.jl
@@ -10,6 +10,13 @@ module QTLLST2DL1
     
     export qtllst2dl1
 
+    """
+    qtllst2dl1(signal, lmax)
+
+    Find a quadtree level-list up to depth `lmax` that minimizes a DST-based cost.
+
+    Returns `(liste, cost)` where `liste` encodes the partition and `cost` is the objective value.
+    """
     function qtllst2dl1(signal::AbstractVecOrMat, lmax::Int)
         liste, cost = recurs_testdst2d(signal, 0, lmax)
         return liste, cost
diff --git a/src/solvelaplace.jl b/src/solvelaplace.jl
index 4138eed..7c4c1d8 100644
--- a/src/solvelaplace.jl
+++ b/src/solvelaplace.jl
@@ -33,6 +33,17 @@ module SOLVELAPLACE
     end
 
 
+    """
+    solvelaplace(image)
+
+    Solve the Laplace (Poisson) equation on the interior of `image` using DST/IDST.
+
+    Arguments
+    - `image`: m×n array with boundary values; the routine solves for the interior (returns array sized m×n with boundary preserved).
+
+    Returns
+    - array with the computed solution where interior values have been solved and boundaries are unchanged.
+    """
     function solvelaplace(image::AbstractVecOrMat)
         if (isempty(image) == true)
             error("IM has to be a 2D array")
@@ -127,7 +138,13 @@ module SOLVELAPLACE
         i1 = idst(hi1)
         i2 = idst(h2m' + h4m')'
         lapl = lapl + i1 + i2
-        #lapl = lapl + idst(h1m+h3m) + idst(h2m'+h4m')';
-        return lapl
+        # Place interior solution back into full image and preserve boundaries
+        out = zeros(m, n)
+        out[2:end-1, 2:end-1] .= lapl
+        out[1, :] .= image[1, :]
+        out[end, :] .= image[end, :]
+        out[:, 1] .= image[:, 1]
+        out[:, end] .= image[:, end]
+        return out
     end
 end;
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index aa50c07..333c256 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,7 +21,63 @@ has_phlct_restore = isdefined(PolyHarmonicTrigTransforms, :phlct_restore);
 has_phlct2d = isdefined(PolyHarmonicTrigTransforms, :phlct2d);
 @test has_phlct2d = true;
 
-has_laplace = isdefined(PolyHarmonicTrigTransforms, :laplace);
-@test has_laplace::Bool = true;
+# DST/IDST round-trip
+x = randn(16)
+y = dst(x)
+x_rec = idst(y)
+@test isapprox(x, x_rec; atol=1e-10, rtol=1e-8)
 
-print("tests passed");
\ No newline at end of file
+# PHLCT forward/backward round-trip
+n = 8
+img = randn(n, n)
+coeffs = phlct_forward(img, n)
+img_rec = phlct_backward(coeffs, n)
+@test isapprox(img, img_rec; atol=1e-10, rtol=1e-8)
+
+# In-place `_!` variant tests
+# DST!/IDST!
+x2 = randn(16)
+xp = copy(x2)
+yp = dst(copy(x2))
+dst!(xp)
+@test isapprox(xp, yp; atol=1e-10, rtol=1e-8)
+
+id_in = dst(randn(16))
+idp = copy(id_in)
+idst!(idp)
+@test isapprox(idst(id_in), idp; atol=1e-10, rtol=1e-8)
+
+# PHLCT `_!` variants
+img2 = randn(n, n)
+coeffs_expected = phlct_forward(img2, n)
+img2p = copy(img2)
+phlct_forward!(img2p, n)
+@test isapprox(img2p, coeffs_expected; atol=1e-10, rtol=1e-8)
+
+coeffs2 = phlct_forward(randn(n, n), n)
+coeffs2p = copy(coeffs2)
+phlct_backward!(coeffs2p, n)
+@test isapprox(phlct_backward(coeffs2, n), coeffs2p; atol=1e-10, rtol=1e-8)
+
+# LLST 2D `_!` variants
+data = randn(12, 12)
+res_ll = llst2d(copy(data))
+data_p = copy(data)
+llst2d!(data_p)
+@test isapprox(data_p, res_ll; atol=1e-10, rtol=1e-8)
+
+res_ill = illst2d(copy(data))
+data_p2 = copy(data)
+illst2d!(data_p2)
+@test isapprox(data_p2, res_ill; atol=1e-10, rtol=1e-8)
+
+# Basic solvelaplace smoke test if available
+if isdefined(PolyHarmonicTrigTransforms, :solvelaplace)
+	im = zeros(6, 6)
+	im[1, :] .= 1.0
+	im[end, :] .= 1.0
+	res = solvelaplace(im)
+	@test size(res) == size(im)
+end
+
+println("All tests passed.")
\ No newline at end of file