diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 59d999f..cb0cb6f 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -43,7 +43,10 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@v1
         with:
-          version: '1'
+          version: '1.9'
+      - uses: julia-actions/cache@v1
+        with:
+          cache-registries: "true"
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-docdeploy@v1
         env:
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index c8b80bb..17b438e 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.8.0"
+julia_version = "1.8.5"
 manifest_format = "2.0"
-project_hash = "a66165f6da81458613938be0d3cf4835e0b4522b"
+project_hash = "c314d07c33c737c2dea5d86662f8061ff527d5f7"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
@@ -18,15 +18,15 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
 [[deps.DocStringExtensions]]
 deps = ["LibGit2"]
-git-tree-sha1 = "5158c2b41018c5f7eb1470d558127ac274eca0c9"
+git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
 uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.9.1"
+version = "0.9.3"
 
 [[deps.Documenter]]
 deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "6030186b00a38e9d0434518627426570aac2ef95"
+git-tree-sha1 = "58fea7c536acd71f3eef6be3b21c0df5f3df88fd"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.27.23"
+version = "0.27.24"
 
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
@@ -46,8 +46,8 @@ version = "0.21.3"
 
 [[deps.JuliaUtils]]
 path = ".."
-uuid = "5394660d-8fe9-41aa-abe5-b2992337cf4a"
-version = "0.1.0"
+uuid = "e738e926-3d3f-4c6c-9019-dd7448bc78fc"
+version = "1.0.0"
 
 [[deps.LibGit2]]
 deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
@@ -68,10 +68,16 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.2.0"
 
 [[deps.Parsers]]
-deps = ["Dates"]
-git-tree-sha1 = "3d5bf43e3e8b412656404ed9466f1dcbf7c50269"
+deps = ["Dates", "SnoopPrecompile"]
+git-tree-sha1 = "478ac6c952fddd4399e71d4779797c538d0ff2bf"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.4.0"
+version = "2.5.8"
+
+[[deps.Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.3.0"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -92,9 +98,20 @@ version = "0.7.0"
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[deps.SnoopPrecompile]]
+deps = ["Preferences"]
+git-tree-sha1 = "e760a70afdcd461cf01a575947738d359234665c"
+uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c"
+version = "1.0.3"
+
 [[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
+[[deps.TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+version = "1.0.0"
+
 [[deps.Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/docs/make.jl b/docs/make.jl
index 4388324..90c5b90 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -15,7 +15,26 @@ makedocs(;
         assets=String[],
     ),
     pages=[
-        "Home" => "index.md",
+    "Home" => "index.md",
+	"Core utility functions" => [ 
+        "I/O"=> "core/io.md",
+        "Space"=> "core/space.md",
+        "Periodic"=> "core/periodic.md",
+        "Physics"=> "core/physics.md",
+        "Defects"=> "core/defects.md",
+        "Spectrum"=> "core/spectrum.md",
+        "Interpolations"=> "core/interpolations.md"
+        ],
+	"Ploting Utilities" => [
+        "Makie"    => "ploting/makie.md",
+        "PyPlot"   => "ploting/pyplot.md"
+        ],
+	"Macros" =>  "macros/macros.md",
+    "Library" => [
+        "Public" => "library/public.md", 
+        "Internals" => "library/internals.md",
+        "Glossary" => "library/glossary.md"
+        ]
     ],
 )
 
diff --git a/docs/src/core/defects.md b/docs/src/core/defects.md
new file mode 100644
index 0000000..2179b76
--- /dev/null
+++ b/docs/src/core/defects.md
@@ -0,0 +1,4 @@
+```@autodocs
+ Modules = [JuliaUtils]
+ Pages   = ["defects.jl"]
+```
diff --git a/docs/src/core/interpolations.md b/docs/src/core/interpolations.md
new file mode 100644
index 0000000..888d219
--- /dev/null
+++ b/docs/src/core/interpolations.md
@@ -0,0 +1,4 @@
+```@autodocs
+ Modules = [JuliaUtils]
+ Pages   = ["interpolations.jl"]
+```
diff --git a/docs/src/core/io.md b/docs/src/core/io.md
new file mode 100644
index 0000000..2cd1eb6
--- /dev/null
+++ b/docs/src/core/io.md
@@ -0,0 +1,4 @@
+```@autodocs
+ Modules = [JuliaUtils]
+ Pages   = ["io.jl"]
+```
diff --git a/docs/src/core/periodic.md b/docs/src/core/periodic.md
new file mode 100644
index 0000000..ca4373f
--- /dev/null
+++ b/docs/src/core/periodic.md
@@ -0,0 +1,4 @@
+```@autodocs
+ Modules = [JuliaUtils]
+ Pages   = ["periodic.jl"]
+```
diff --git a/docs/src/core/physics.md b/docs/src/core/physics.md
new file mode 100644
index 0000000..a40692d
--- /dev/null
+++ b/docs/src/core/physics.md
@@ -0,0 +1,4 @@
+```@autodocs
+ Modules = [JuliaUtils]
+ Pages   = ["physics.jl"]
+```
diff --git a/docs/src/core/space.md b/docs/src/core/space.md
new file mode 100644
index 0000000..0d29594
--- /dev/null
+++ b/docs/src/core/space.md
@@ -0,0 +1,4 @@
+```@autodocs
+ Modules = [JuliaUtils]
+ Pages   = ["space.jl"]
+```
diff --git a/docs/src/core/spectrum.md b/docs/src/core/spectrum.md
new file mode 100644
index 0000000..35b429d
--- /dev/null
+++ b/docs/src/core/spectrum.md
@@ -0,0 +1,4 @@
+```@autodocs
+ Modules = [JuliaUtils]
+ Pages   = ["spectrum.jl"]
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index af22d20..6cb9881 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -4,11 +4,20 @@ CurrentModule = JuliaUtils
 
 # JuliaUtils
 
-Documentation for [JuliaUtils](https://github.com/navdeeprana/JuliaUtils.jl).
+Documentation for [JuliaUtils](https://github.com/pankajpopli/JuliaUtils.jl).
+
+## Overview
+
+## Instalation
+
+Install this package by using
+add github_repo_url
+
+
+## Devlopers
+
+This repo is currently being devloped by Navdeep Rana and Pankaj Popli.
+
+## Cite
 
-```@index
-```
 
-```@autodocs
-Modules = [JuliaUtils]
-```
diff --git a/docs/src/library/glossary.md b/docs/src/library/glossary.md
new file mode 100644
index 0000000..69a2afe
--- /dev/null
+++ b/docs/src/library/glossary.md
@@ -0,0 +1,12 @@
+## Modules
+
+```@autodocs
+Modules = [JuliaUtils,JuliaUtils.MakieUtils,JuliaUtils.PyPlotUtils]
+Order  = [:module]
+```
+
+
+## Index of all functions
+```@index
+Modules = [JuliaUtils,JuliaUtils.PyPlotUtils,JuliaUtils.MakieUtils]
+```
diff --git a/docs/src/library/internals.md b/docs/src/library/internals.md
new file mode 100644
index 0000000..090524c
--- /dev/null
+++ b/docs/src/library/internals.md
@@ -0,0 +1,3 @@
+## Mention any internal details here
+
+I have created this file in spirit of other published documentations. Should remove if not needed.
diff --git a/docs/src/library/public.md b/docs/src/library/public.md
new file mode 100644
index 0000000..134bb0f
--- /dev/null
+++ b/docs/src/library/public.md
@@ -0,0 +1,3 @@
+## Any public details about the package
+
+I have created this file in spirit of other published documentations. Should remove if not needed.
diff --git a/docs/src/macros/macros.md b/docs/src/macros/macros.md
new file mode 100644
index 0000000..ec786cf
--- /dev/null
+++ b/docs/src/macros/macros.md
@@ -0,0 +1,7 @@
+## Custom macros definitions
+
+```@autodocs
+Modules = [JuliaUtils,JuliaUtils.MakieUtils,JuliaUtils.PyPlotUtils]
+Order = [:macro]
+Pages   = ["macros.jl"]
+```
diff --git a/docs/src/ploting/makie.md b/docs/src/ploting/makie.md
new file mode 100644
index 0000000..fad3a6c
--- /dev/null
+++ b/docs/src/ploting/makie.md
@@ -0,0 +1,5 @@
+```@autodocs
+ Modules = [JuliaUtils.MakieUtils]
+Pages   = ["makie.jl"]
+Order = [:constant,:type,:function,:macro]
+```
diff --git a/docs/src/ploting/pyplot.md b/docs/src/ploting/pyplot.md
new file mode 100644
index 0000000..b0dea5d
--- /dev/null
+++ b/docs/src/ploting/pyplot.md
@@ -0,0 +1,5 @@
+```@autodocs
+ Modules = [JuliaUtils.PyPlotUtils]
+ Pages   = ["pyplot.jl"]
+Order = [:constant,:type,:function,:macro]
+```
diff --git a/src/JuliaUtils.jl b/src/JuliaUtils.jl
index 209d74c..698ca8f 100644
--- a/src/JuliaUtils.jl
+++ b/src/JuliaUtils.jl
@@ -1,54 +1,72 @@
 """
-Main module for `JuliaUtils.jl` -- my julia language utility functions.
-
+Main module for `Julia Useful functions`
 # Exports
 $(EXPORTS)
 """
 module JuliaUtils
+   #core exported functions
+   export read_data,
+      read_xydata,
+      write_data,
+      write_xydata,
+      linspace,
+      logspace,
+      logspacedPick,
+      logspacedPick!,
+      cartesian_mesh,
+      polar_mesh,
+      periodic_distance,
+      aperiodic,
+      defect_vector_field,
+      find_defects,
+      chargeAt,
+      nearest_neighbour_distance,
+      radial_distribution_function,
+      Interp2D,
+      order_parameter,
+      @repeat
 
-export read_data,
-    read_xydata,
-    write_data,
-    write_xydata,
-    linspace,
-    logspaced,
-    cartesian_mesh,
-    polar_mesh,
-    periodic_distance,
-    aperiodic,
-    defect_vector_field,
-    find_defects,
-    charge,
-    nearest_neighbour_distance,
-    radial_distribution_function,
-    @repeat
-
-using DocStringExtensions
-
-include("io.jl")
-include("space.jl")
-include("defects.jl")
-include("physics.jl")
-include("macros.jl")
-
-module PyPlotUtils
+   using DocStringExtensions
+   include("io.jl")
+   include("space.jl")
+   include("defects.jl")
+   include("physics.jl")
+   include("macros.jl")
+   include("interpolations.jl")
+   include("periodic.jl")
+   include("spectrum.jl")
+   include("pyplot.jl")
+   include("makie.jl")
 
-export colors,
-    figax,
-    noticks!,
-    fontsize!,
-    formatter!
+   # include all files recursively from the `src` folder
+   # for (root, dirs, files) in walkdir("src")
+   #    for file in files
+   #       # check if the file is a Julia file
+   #       if endswith(file, ".jl")
+   #          # construct the path to the file
+   #          path = joinpath(root, file)
+   #          # include the file in the module
+   #          Base.include(path)
+   #       end
+   #    end
+   # end
 
-include("pyplot.jl")
 
-end
-
-module MakieUtils
 
-export figax
 
-include("makie.jl")
-
-end
+   # module PyPlotUtils
+   # export colors,
+   #    figax,
+   #    noticks!,
+   #    fontsize!,
+   #    formatter!
+   #    include("pyplot.jl")
+   # end
+   
+   # module MakieUtils
+   # using DocStringExtensions
+   # export figax
+   #    include("makie.jl")
+   # end
 
 end
diff --git a/src/defects.jl b/src/defects.jl
index 0dde920..c20725e 100644
--- a/src/defects.jl
+++ b/src/defects.jl
@@ -1,102 +1,168 @@
 """
-Generate a topological defect vector field in two dimensions.
+Generate a single topological defect vector field in two dimensions.
+$(TYPEDSIGNATURES) 
+Args:\\
+[Req] q::Integer   (value of the charge)\\
+[Req] lx::Real     (linear size of the box. The square box is used.)\\
+[Req] m::Integer   (value of topological charge.)\\
+[Req] nx::Integer  (number of grid points in the box, with ny=nx)\\
+[Opt] θ0 = 0       (constant overall global phase, default is zero.)\\
+[Opt] kwargs...    (extra keyword arguments for cartesian_mesh() function.)
+Note:\\
+functionn generates its own box and cartesian mesh. Use other method to provide mesh coordinates explictly.
+@Examples
+a = 2
+b =4 
+a+b
+"""
+function defect_vector_field(q::Integer, lx::Real, nx::Integer; θ0 = 0, kwargs...)
+   x, y = cartesian_mesh(nx, lx; kwargs...)
+   θ = @. atan(y, x)
+   u = @. cos(q * θ + θ0)
+   v = @. sin(q * θ + θ0)
+   return u, v
+end
+"""
+Generate vector field due to many topological defectsin two dimensions.
+$(TYPEDSIGNATURES) 
+Args:\\
+[Req] x::AbstractArray{U} (x coordinates, in XY Form, xy' Form or mesh Form)\\
+[Req] y::AbstractArray{U} (x coordinates, in XY Form, xy' Form or mesh Form)\\
+[Req] Q::AbstractVector{<:Real} (Array of charge values,[q1,q2,...])\\
+[Req] postQ::AbstractVector{T}  (Array of tuples for coordinates of charges, [(ax1,ay1),(ax2,ay2),...])\\
+[Opt] θ0 = 0 (global phase for the vector field)
+"""
+function defect_vector_field(
+   x::AbstractArray{U},
+   y::AbstractArray{U},
+   Q::AbstractVector{<:Real},
+   postQ::AbstractVector{T};
+   θ0 = 0) where {T<:Tuple{Real,Real},U<:Real}
+
+   @assert (typeof(x)<:AbstractVector && y'==adjoint(x) ? false : true)  "expected (x,y <: AbstractVector{Real} in XY form) or (x,y <:AbstractArray{Real,2} in mesh form) or (y<:Adjoint of x)"
+   @assert (size(Q)[1] == size(postQ)[1]) "Number of charges does not have equivalent postion array."
 
-$(SIGNATURES)
-"""
-function defect_vector_field(m::Integer, lx, nx; θ0 = 0, kwargs...)
-    x, y = cartesian_mesh(nx, lx; kwargs...)
-    θ = @. atan(y, x)
-    u = @. cos(m * θ + θ0)
-    v = @. sin(m * θ + θ0)
-    return u, v
+   θ = fill!(similar(x.*y),0.0)
+   @inbounds for (q,(ax,ay)) in zip(Q,postQ)
+     @. θ +=  q*atan(y-ay, x-ax)
+   end
+   u = @. cos(θ + θ0)
+   v = @. sin(θ + θ0)
+   return u,v
 end
 
 @inline function _get_next(i::T, N::T) where {T<:Integer}
-    i == N ? 1 : i + 1
+   i == N ? 1 : i + 1
 end
 
 @inline function _area(tn::T, t0::T) where {T<:Real}
-    dt = tn - t0
-    if (dt > π)
-        dt = dt - 2π
-    end
-    if (dt < -π)
-        dt = dt + 2π
-    end
-    return dt
+   dt = tn - t0
+   if (dt > π)
+      dt = dt - 2π
+   end
+   if (dt < -π)
+      dt = dt + 2π
+   end
+   return dt
 end
 
+"""
+$(METHODLIST)
+"""
+function find_defects end
 """
 Count all topological zeros for a 2D orientation field.
-
-$(SIGNATURES)
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] θ::AbstractArray{T,2}; (Array of θ values, make sure they are column ordered)\\
+[Opt] θ0::Real = 0.8 * (2π) (optional threshold value above which a charge is counted. Doesn't affect much.)\\
+[Opt] mx::Integer = 0 (scan in a small window of the box, anchored at (0,0) and streched upto (nx-mx,ny-my))\\
+[Opt] my::Integer = mx (scan in a small window of the box, anchored at (0,0) and streched upto (nx-mx,ny-my))
 """
-function find_defects(θ::AbstractArray{T,2}; θ0 = 0.8 * (2π), mx = 0, my = mx) where {T<:Real}
-    nx, ny = size(θ)
-    x, y, c = Int[], Int[], Int[]
+function find_defects(θ::AbstractArray{T,2}; θ0::Real = 0.8 * (2π), mx::Integer = 0, my::Integer = mx) where {T<:Real}
+   nx, ny = size(θ)
+   x, y, c = Int[], Int[], Int[]
 
-    @inbounds for j in 1:ny-my
-        jnext = _get_next(j, ny)
-        @inbounds for i in 1:nx-mx
-            inext = _get_next(i, nx)
-            t1 = θ[i, j]
-            t2 = θ[inext, j]
-            t3 = θ[inext, jnext]
-            t4 = θ[i, jnext]
+   @inbounds for j in 1:ny-my
+      jnext = _get_next(j, ny)
+      @inbounds for i in 1:nx-mx
+         inext = _get_next(i, nx)
+         t1 = θ[i, j]
+         t2 = θ[inext, j]
+         t3 = θ[inext, jnext]
+         t4 = θ[i, jnext]
 
-            Δθ = _area(t2, t1) + _area(t3, t2) + _area(t4, t3) + _area(t1, t4)
+         Δθ = _area(t2, t1) + _area(t3, t2) + _area(t4, t3) + _area(t1, t4)
 
-            if abs(Δθ) > θ0
-                push!(x, i)
-                push!(y, j)
-                push!(c, round(Int, Δθ / (2π)))
-            end
-        end
-    end
-    return x, y, c
+         if abs(Δθ) > θ0
+               push!(x, i)
+               push!(y, j)
+               push!(c, round(Int, Δθ / (2π)))
+         end
+      end
+   end
+   return x, y, c
 end
-
+"""
+Count all topological zeros for a 2D orientation field given cos and sin values.
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] u::AbstractArray{T,2}; (Array of cos(θ) values, make sure they are column ordered)\\
+[Req] u::AbstractArray{T,2}; (Array of sin(θ) values, make sure they are column ordered)\\
+[Opt] kwargs... (Same as find_defects() function)\\
+"""
 function find_defects(u::AbstractArray{T,2}, v::AbstractArray{T,2}; kwargs...) where {T<:Real}
-    θ = @. atan(v, u)
-    return find_defects(θ; kwargs...)
+   θ = @. atan(v, u)
+   return find_defects(θ; kwargs...)
 end
-
+"""
+Count all topological zeros for a 2D orientation field given complex field values.
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] ϕ::AbstractArray{Complex{T},2} (Array of cos(θ) + ι sin(θ), column ordered)\\
+[Opt] kwargs... (Same as find_defects() function)\\
+"""
 function find_defects(ϕ::AbstractArray{Complex{T},2}; kwargs...) where {T<:Real}
-    θ = @. angle(ϕ)
-    return find_defects(θ; kwargs...)
+   θ = @. angle(ϕ)
+   return find_defects(θ; kwargs...)
 end
 
 """
 Compute topological charge at points (x,y) for the orientation field θ
-
-$(SIGNATURES)
-"""
-function charge(x, y, θ::AbstractArray{T,2}; Δ = 1) where {T<:Real}
-    nx, ny = size(θ)
-    c = zeros(Float64, length(x))
-    for (n, (xn, yn)) in enumerate(zip(x, y))
-        for Δy in -Δ:+Δ
-            for Δx in -Δ:+Δ
-                i, j = mod(xn + Δx - 1, nx) + 1, mod(yn + Δy - 1, ny) + 1
-                inext = _get_next(i, nx)
-                jnext = _get_next(j, ny)
-                t1 = θ[i, j]
-                t2 = θ[inext, j]
-                t3 = θ[inext, jnext]
-                t4 = θ[i, jnext]
-                c[n] += (_area(t2, t1) + _area(t3, t2) + _area(t4, t3) + _area(t1, t4)) / (2π)
-            end
-        end
-    end
-    return c
-end
-
-function charge(x, y, u::AbstractArray{T,2}, v::AbstractArray{T,2}; kwargs...) where {T<:Real}
-    θ = @. atan(v, u)
-    return charge(x, y, θ; kwargs...)
+$(TYPEDSIGNATURES)
+"""
+function chargeAt(x, y, θ::AbstractArray{T,2}; Δ = 1) where {T<:Real}
+   nx, ny = size(θ)
+   c = zeros(Float64, length(x))
+   for (n, (xn, yn)) in enumerate(zip(x, y))
+      for Δy in -Δ:+Δ
+         for Δx in -Δ:+Δ
+               i, j = mod(xn + Δx - 1, nx) + 1, mod(yn + Δy - 1, ny) + 1
+               inext = _get_next(i, nx)
+               jnext = _get_next(j, ny)
+               t1 = θ[i, j]
+               t2 = θ[inext, j]
+               t3 = θ[inext, jnext]
+               t4 = θ[i, jnext]
+               c[n] += (_area(t2, t1) + _area(t3, t2) + _area(t4, t3) + _area(t1, t4)) / (2π)
+         end
+      end
+   end
+   return c
 end
-
-function charge(x, y, ϕ::AbstractArray{Complex{T},2}; kwargs...) where {T<:Real}
-    θ = @. angle(ϕ)
-    return charge(x, y, θ; kwargs...)
+"""
+Compute topological charge at points (x,y) for the vector field cos(θ),sin(θ)
+$(TYPEDSIGNATURES)
+"""
+function chargeAt(x, y, u::AbstractArray{T,2}, v::AbstractArray{T,2}; kwargs...) where {T<:Real}
+   θ = @. atan(v, u)
+   return chargeAt(x, y, θ; kwargs...)
 end
+"""
+Compute topological charge at points (x,y) for the complex field cos(θ)+ ι sin(θ)
+$(TYPEDSIGNATURES)
+"""
+function chargeAt(x, y, ϕ::AbstractArray{Complex{T},2}; kwargs...) where {T<:Real}
+   θ = @. angle(ϕ)
+   return chargeAt(x, y, θ; kwargs...)
+end
\ No newline at end of file
diff --git a/src/interpolations.jl b/src/interpolations.jl
new file mode 100644
index 0000000..05f9d57
--- /dev/null
+++ b/src/interpolations.jl
@@ -0,0 +1,26 @@
+using Interpolations
+"""
+Interpolate a 2D data 
+$(TYPEDSIGNATURES)
+- [Req] data::AbstractArray{T,1} (input data vector, row ordered (0,0) at bottom left corner)
+- [Req] order::T=1.3   (order of interpolation)
+Note:\\
+data is strutred with (row+nx*col) rule starting from bottom left corner.
+"""
+function Interp2D(data::AbstractArray{T,1}; order::T=1.3) where T<:Real
+    
+   IC = CubicSplineInterpolation((axes(data,1), axes(data,2)), data)
+
+   finerx = LinRange(firstindex(data,1), lastindex(data,1), size(data,1) * order)
+   finery = LinRange(firstindex(data,2), lastindex(data,2), size(data,2) * order)
+   nx = length(finerx)
+   ny = length(finery)
+
+   data_interp = Array{Float64}(undef,nx,ny)
+   for i ∈ 1:nx, j ∈ 1:ny
+      data_interp[i,j] = IC(finerx[i],finery[j])
+   end
+
+   return finerx, finery, data_interp
+
+end
diff --git a/src/io.jl b/src/io.jl
index 89fcdc4..19f26c0 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -1,42 +1,120 @@
 using CSV
 using DataFrames
 using HDF5
+using DelimitedFiles
 
+"""
+$(TYPEDSIGNATURES)
+!!! note "Efficiency"
+    It is noted that CSV is faster compared to Readdlm package for large files.
+    The loadtxt function defined below uses readdlm.
+    Check the timings using @time macro for your needs.
+    see ["stackoverflow"](https://stackoverflow.com/questions/26810171/is-there-an-equivalent-or-close-to-numpy-loadtxt-for-julia)
+
+!!! warning "Not tested"
+    Read-write using CSV and HDF5 are not tested fully. Use [`loadtxt`](@ref) function for the time being.\\
+
+Read from HDF5 file.
+Args:\\
+[Req] fname; (filename)
+[Opt] slice = (:, :) (Read HDF5 documentation)\\
+[Opt] dsets = ? (Read HDF5 documentation)
+"""
 function read_data(fname; slice = (:, :), dsets = ["phi1", "phi2"])
-    return [h5read(fname, dset, slice) for dset in dsets]
+   return [h5read(fname, dset, slice) for dset in dsets]
 end
+"""
 
+$(TYPEDSIGNATURES)
+Write to HDF5 file.
+Args:\\
+[Req] fname; (filename)\\
+[Opt] mode = "cw", (Read HDF5 documentation)\\
+[Opt] overwrite = true, (Read HDF5 documentation)\\
+[Opt] kwargs...   (Read HDF5 documentation)\\
+"""
 function write_data(fname; mode = "cw", overwrite = true, kwargs...)
-    h5open(fname, mode) do fid
-        for (k, v) in pairs(kwargs)
-            ds = string(k)
-            if haskey(fid, ds) & overwrite
-                delete_object(fid, ds)
-            end
-            fid[ds] = v
-        end
-    end
+   h5open(fname, mode) do fid
+      for (k, v) in pairs(kwargs)
+         ds = string(k)
+         if haskey(fid, ds) & overwrite
+               delete_object(fid, ds)
+         end
+         fid[ds] = v
+      end
+   end
 end
-
+"""
+$(TYPEDSIGNATURES)
+Read and XY format data file usig CSV package.
+Args:\\
+[Req] fname; (filename)\\
+[Opt] delim = " ",   (read documentation for CSV.read)
+[Opt] header = false,   (read documentation for CSV.read)
+[Opt] fcol = :Column1,  (read documentation for CSV.read)
+[Opt] ffun = (x -> true),  (read documentation for CSV.read)
+[Opt] kwargs...   (read documentation for CSV.read)
+"""
 function read_xydata(fname; delim = " ", header = false, fcol = :Column1, ffun = (x -> true), kwargs...)
-    df = CSV.read(
-        fname,
-        DataFrame;
-        delim = delim,
-        ignorerepeated = true,
-        header = header,
-        ntasks = 1,
-        kwargs...
-    )
-    return eachcol(filter(fcol => ffun, df))
+   df = CSV.read(
+      fname,
+      DataFrame;
+      delim = delim,
+      ignorerepeated = true,
+      header = header,
+      ntasks = 1,
+      kwargs...
+   )
+   return eachcol(filter(fcol => ffun, df))
 end
-
+"""
+$(TYPEDSIGNATURES)
+Wrtie XY format data to a file using CSV package.\\
+Args:\\
+[Req] fname (filename)\\
+[Opt] data (data to be written)
+[Opt] kwargs...   (read documentation for CSV.read)
+"""
 function write_xydata(fname, data; kwargs...)
-    return CSV.write(
-        fname,
-        DataFrame(data, :auto),
-        delim = " ",
-        writeheader = false,
-        kwargs...
-    )
+   return CSV.write(
+      fname,
+      DataFrame(data, :auto),
+      delim = " ",
+      writeheader = false,
+      kwargs...
+   )
+end
+"""
+$(TYPEDSIGNATURES)
+Read XY format data file using Readdlm package.\\
+Args:\\
+[Req] fname (filename)\\
+[Opt] skiprow (# of header to be skiped)\\
+[Opt] delim::AbstractChar = '\t',   (seperator)\\
+[Opt] eol::AbstractChar='\n', (end of line character)\\
+[Opt] kwargs...   (read documentation forDelimitedFiles.readdlm)
+"""
+function loadtxt(fname; skiprows=0, delim::AbstractChar = '\t',eol::AbstractChar='\n', kwargs...)
+   return readdlm(
+      fname,
+      delim,
+      eol,
+      skipstart=skiprows,
+      skipblanks = true,
+      kwargs...)
+end
+"""
+$(TYPEDSIGNATURES)
+Write XY format data file using Readdlm package.\\
+Args:\\
+ - `fname` (filename)\\
+ - `data` (data to be written, Vector/Matrix or an iterable collection of iterable rows)\\
+ - `mode = "w"` (mode of opening file, "w" = opwn for write)
+ - `delim::AbstractChar = '\t'`   (seperator)\\
+ - `kwargs...`   (read documentation for DelimitedFiles.writedlm)\\
+"""
+function savetxt(fname, data; mode="w", delim::AbstractChar = '\t', kwargs...)
+   open(fname, mode) do fid
+      writedlm(fid, data,delim,kwargs...)
+  end
 end
diff --git a/src/macros.jl b/src/macros.jl
index b77e368..fe92c10 100644
--- a/src/macros.jl
+++ b/src/macros.jl
@@ -1,18 +1,22 @@
 # Copied from https://jkrumbiegel.com/pages/2021-06-07-macros-for-beginners/
+"""
+$(TYPEDSIGNATURES)
+@macro repeat(expr, sizes...)
+"""
 macro repeat(exp, sizes...)
-    iterator_expressions = map(sizes) do s
-        Expr(
-            :(=),
-            :_,
-            quote
-                1:$(esc(s))
-            end
-        )
-    end
+   iterator_expressions = map(sizes) do s
+       Expr(
+           :(=),
+           :_,
+           quote
+               1:$(esc(s))
+           end
+       )
+   end
 
-    Expr(
-        :comprehension,
-        esc(exp),
-        iterator_expressions...
-    )
-end
+   Expr(
+       :comprehension,
+       esc(exp),
+       iterator_expressions...
+   )
+end
\ No newline at end of file
diff --git a/src/makie.jl b/src/makie.jl
index 96d68ed..bb8bb84 100644
--- a/src/makie.jl
+++ b/src/makie.jl
@@ -1,12 +1,27 @@
-using Makie: Figure, Axis, AxisAspect
+"""
+Makie module for JuliaUseful Ploting functions Utilities`
+# Exports
+$(EXPORTS)
+"""
+module MakieUtils
 
-function figax(; nx = 1, ny = 1, h = 3, a = 1.6, s = 100, kwargs...)
-    res = (round(Int, a * s * h * nx), round(Int, s * h * ny))
-    fig = Figure(resolution = res)
-    ax = [Axis(fig[j, i]; aspect = AxisAspect(a), kwargs...) for i in 1:nx, j in 1:ny]
-    if nx * ny == 1
-        return fig, ax[1]
-    else
-        return fig, ax
-    end
-end
+   using Makie: Figure, Axis, AxisAspect
+   using DocStringExtensions
+
+   export figax
+
+   """
+   $(TYPEDSIGNATURES)
+   Create Figure, Axis plots for input number of row and columns.
+   """
+   function figax(; nx::Int = 1, ny::Int = 1, h = 3, a = 1.6, s = 100, kwargs...)
+      res = (round(Int, a * s * h * nx), round(Int, s * h * ny))
+      fig = Figure(resolution = res)
+      ax = [Axis(fig[j, i]; aspect = AxisAspect(a), kwargs...) for i in 1:nx, j in 1:ny]
+      if nx * ny == 1
+         return fig, ax[1]
+      else
+         return fig, ax
+      end
+   end
+end #module
\ No newline at end of file
diff --git a/src/periodic.jl b/src/periodic.jl
new file mode 100644
index 0000000..905bd7d
--- /dev/null
+++ b/src/periodic.jl
@@ -0,0 +1,124 @@
+"""
+Wraps angle θ in [-π,π) or [0,2π) domain.\\
+$(TYPEDSIGNATURES)
+[Req] θ::Real (input angle value)
+[Opt] domain=π (domain type, π => [-π,π) and 2π =  [0,2π) domain)
+"""
+function wrap_angle(θ::Real;domain=π)
+   @assert domain==π || domain==2π "'domain' can take values of π or 2π, default is π"
+   θ = mod(θ + domain,2π)
+   if (θ < 0)
+     θ += 2π
+   end
+   return θ - domain
+end
+"""
+Wraps angle difference (θ2 - θ1) in [-π,π) or [0,2π) domain.\\
+benchmark with _area(θ2,θ1) before using this function.
+$(TYPEDSIGNATURES)
+[Req] θ1::Real (input angle1 value)
+[Req] θ2::Real (input angle2 value)
+[Opt] domain=π (domain type, π => [-π,π) and 2π =  [0,2π) domain)
+"""
+function wrap_angle(θ2::T,θ1::T;domain=π) where {T<:Real}
+   diff = θ2-θ1
+   return anglewrap(diff,domain=domain)
+end
+# function anglewrap(angle::Real)
+#    angle = mod(angle + π,2π);
+#    if (angle < 0)
+#       angle += 2π;
+#    end
+#    return angle - π;
+# end
+
+#####################
+"""
+find periodic distance provided absolute distance abs(xi - xj).
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] x::T,  (absolute distance)
+[Req] l::T,  (linear box size)
+[Req] hl::T  (half length of box size)
+"""
+@inline function periodic_distance(x::Real, l::Real, hl::Real)
+   @assert !(x < 0.0) "provide absolute distance between two coordinates abs(xj - xi)"
+   @. ifelse(abs(x) > hl, l - x, x)
+end
+"""
+find periodic distance provided absolute distance abs(xi - xj).
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] x::T,  (absolute distance)\\
+[Req] l::T,  (linear box size)\\
+"""
+@inline function periodic_distance(x::Real, l::Real)
+   @assert !(x<0.0) "provide absolute distance between two coordinates abs(xj - xi)"
+   return periodic_distance(x, l, 0.5 * l)
+end
+"""
+Wrap particle coordinate in [-hl,hl) open range with origin at center.
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] x::T,  (absolute distance)\\
+[Req] l::T,  (linear box size)\\
+[Req] hl::T  (half length of box size)
+"""
+@inline function wrap_cords(x::Real, l::Real, hl::Real)
+   if (abs(x)>hl)
+      if (x >= hl)
+         x = x - l
+      end
+      if (x < -hl)
+         x = x + l
+      end
+   end
+   return x
+end
+"""
+Wrap particle coordinate in [-hl,hl) open range with origin at center.
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] x::T,  (absolute distance)\\
+[Req] l::T,  (linear box size)
+"""
+@inline function wrap_cords(x::Real, l::Real)
+   return wrap_cords(x,l,0.5*l)
+end
+# """
+# Wrap particle coordinate in [-hl,hl] closed range with origin at center.
+# $(TYPEDSIGNATURES)
+# Args:\\
+# [Req] x::T,  (absolute distance)\\
+# [Req] l::T,  (linear box size)
+# """
+# function periodic_cord(x::T, l::T, hl::T) where {T<:Real}
+#    if abs(x) > hl
+#       x -= l * round(x/l)
+#    end
+#    return x
+# end
+#############################################################
+"""
+Convert periodic data to an aperiodic one.
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] x::Vector{T};  (input data vector)\\
+[Opt] period::Real = 2π
+"""
+function aperiodic(x::AbstractVector{T}; period::Real = 2π) where {T<:Real}
+   y = zero(x)
+   x0, Δ, Δp = 0, 0, 0.9 * period
+   for (i, xi) in enumerate(x)
+      Δx = xi - x0
+      if (abs(Δx) >= Δp)
+         Δ = Δ - sign(Δx) * period
+      end
+      y[i] = xi + Δ
+      x0 = xi
+   end
+   return y
+end
+#####################
+
+
diff --git a/src/physics.jl b/src/physics.jl
index a58447a..bdd3589 100644
--- a/src/physics.jl
+++ b/src/physics.jl
@@ -1,57 +1,73 @@
 """
 Compute nearest neighbour distance for a given set of points in a 2D periodic box.
-
-$(SIGNATURES)
+Useful for finding topological nearest neighbours.
+$(TYPEDSIGNATURES)
 """
-function nearest_neighbour_distance(x::Array{T}, y::Array{T}, L::T) where {T<:Real}
-    N = length(x)
-    d = zero(x)
-    hL = 0.5 * L
-    @inbounds for i in 1:N
-        nnd = typemax(T)
-        @inbounds for j in 1:N
-            if (i == j)
-                continue
-            end
-            dx = periodic_distance(abs(x[i] - x[j]), L, hL)
-            dy = periodic_distance(abs(y[i] - y[j]), L, hL)
-            dr = dx^2 + dy^2
-            if (dr < nnd)
-                nnd = dr
-            end
-        end
-        d[i] = sqrt(nnd)
-    end
-    return d
+function nearest_neighbour_distance(x::AbstractArray{T}, y::AbstractArray{T}, L::T) where {T<:Real}
+   N = length(x)
+   d = zero(x)
+   hL = 0.5 * L
+   @inbounds for i in 1:N
+      nnd = typemax(T)
+      @inbounds for j in 1:N
+         if (i == j)
+               continue
+         end
+         dx = periodic_distance(abs(x[i] - x[j]), L, hL)
+         dy = periodic_distance(abs(y[i] - y[j]), L, hL)
+         dr = dx^2 + dy^2
+         if (dr < nnd)
+               nnd = dr
+         end
+      end
+      d[i] = sqrt(nnd)
+   end
+   return d
 end
 
 """
 Compute radial distribution function for a given set of points in a 2D periodic box.
+$(TYPEDSIGNATURES)
+"""
+function radial_distribution_function(x::Array{T}, y::Array{T}, L::T; nbin = 256) where {T<:Real}
+   N = length(x)
+   hL = 0.5 * L
+   dbin = L / nbin
+   r = dbin .* linspace(0, nbin, nbin)
+   gr = zeros(nbin)
+   @inbounds for i in 1:N
+      @inbounds for j in i+1:N
+         dx = periodic_distance(abs(x[i] - x[j]), L, hL)
+         dy = periodic_distance(abs(y[i] - y[j]), L, hL)
+
+         dr = sqrt(dx^2 + dy^2)
+         ir = Int(ceil(dr / dbin))
+         gr[ir] = gr[ir] + 1
+      end
+   end
+   rho = (1 / (L^2)) * (N * (N - 1) / 2)
+   @. gr = gr / (rho * 2 * π * r * dbin)
+   return r, gr
+end
 
-$(SIGNATURES)
 """
-function radial_distribution_function(
-    x::Array{T},
-    y::Array{T},
-    L::T;
-    nbin = 256
-) where {T<:Real}
-    N = length(x)
-    hL = 0.5 * L
-    dbin = L / nbin
-    r = dbin .* linspace(0, nbin, nbin)
-    gr = zeros(nbin)
-    @inbounds for i in 1:N
-        @inbounds for j in i+1:N
-            dx = periodic_distance(abs(x[i] - x[j]), L, hL)
-            dy = periodic_distance(abs(y[i] - y[j]), L, hL)
+$(TYPEDSIGNATURES)
+Compute order parameter given cartesian components of the field.\\
+Args:\\
 
-            dr = sqrt(dx^2 + dy^2)
-            ir = Int(ceil(dr / dbin))
-            gr[ir] = gr[ir] + 1
-        end
-    end
-    rho = (1 / (L^2)) * (N * (N - 1) / 2)
-    @. gr = gr / (rho * 2 * π * r * dbin)
-    return r, gr
+   - `u::Vararg{Array{T}}` [cartesian components of the field, accepts arrays of `Float64` or `ComplexF64`]\\
+
+Examples:\\
+
+ - for 3D vector field ` 𝐩 = (u,v,w)`; call the function as `order_parameter(u,v,w)` \\
+ - for tensorial field `𝐐 = [u v; w z]`;  call the function as `order_parameter(u,v,w,z)` \\
+ - Output is a tuple of `(|𝐦|, [⟨u⟩,⟨v⟩,⟨w⟩...])`, where `|𝐦|` is Frobenius norm of input tensor `𝐐` or vector `𝐩`.\\
+"""
+function order_parameter(u::Vararg{Array{T}}) where {T<:Union{Float64,ComplexF64}}
+   op = zeros(T,length(u))
+   for (i, ui) in enumerate(u)
+      op[i] = sum(ui)/length(ui)
+   end
+   return hypot(op...),op
 end
+# https://math.stackexchange.com/questions/1958882/on-the-generalization-of-polar-coordinates-for-n-dimensions
\ No newline at end of file
diff --git a/src/pyplot.jl b/src/pyplot.jl
index 1092238..b97a231 100644
--- a/src/pyplot.jl
+++ b/src/pyplot.jl
@@ -1,42 +1,71 @@
-using IJulia
-using PyCall
-using PyPlot: plt, Figure, matplotlib
+"""
+PyPlot module for JuliaUseful Ploting functions Utilities`
+# Exports
+$(EXPORTS)
+"""
+module PyPlotUtils
+   using IJulia
+   using PyCall
+   using PyPlot: plt, Figure, matplotlib
+   using DocStringExtensions
 
-colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
+   export figax, 
+      noticks!, 
+      fontsize!, 
+      formatter!
+   
+   colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
 
-function IJulia.metadata(x::Figure)
-    w, h = (x.get_figwidth(), x.get_figheight()) .* x.get_dpi()
-    return Dict("image/png" => Dict("width" => w / 2, "height" => h / 2))
-end
+   """
+   $(TYPEDSIGNATURES)
+   """
+   function IJulia.metadata(x::Figure)
+      w, h = (x.get_figwidth(), x.get_figheight()) .* x.get_dpi()
+      return Dict("image/png" => Dict("width" => w / 2, "height" => h / 2))
+   end
 
-function figax(; nx = 1, ny = 1, h = 3, a = 1.6, kwargs...)
-    figsize = (a * h * nx, h * ny)
-    fig, ax = plt.subplots(ny, nx; figsize = figsize, kwargs...)
-    try
-        ax = permutedims(ax, (2, 1))
-    catch e
-        nothing
-    end
-    nx * ny > 1 ? ax = hcat(ax...) : nothing
-    return fig, ax
-end
+   """
+   $(TYPEDSIGNATURES)
+   Create Figure, Axis plots for input number of row and columns.
+   """
+   function figax(; nx = 1, ny = 1, h = 3, a = 1.6, kwargs...)
+      figsize = (a * h * nx, h * ny)
+      fig, ax = plt.subplots(ny, nx; figsize = figsize, kwargs...)
+      try
+         ax = permutedims(ax, (2, 1))
+      catch e
+         nothing
+      end
+      nx * ny > 1 ? ax = hcat(ax...) : nothing
+      return fig, ax
+   end
 
-function noticks!(ax; x = false, y = false)
-    ax.xaxis.set_visible(x)
-    ax.yaxis.set_visible(y)
-    return nothing
-end
+   """
+   $(TYPEDSIGNATURES)
+   """
+   function noticks!(ax; x = false, y = false)
+      ax.xaxis.set_visible(x)
+      ax.yaxis.set_visible(y)
+      return nothing
+   end
 
-function fontsize!(ax, fontsize)
-    for textobj in ax.findobj(match = (x -> x.__module__ == "matplotlib.text"))
-        textobj.set_fontsize(fontsize)
-    end
-    return nothing
-end
+   """
+   $(TYPEDSIGNATURES)
+   """
+   function fontsize!(ax, fontsize)
+      for textobj in ax.findobj(match = (x -> x.__module__ == "matplotlib.text"))
+         textobj.set_fontsize(fontsize)
+      end
+      return nothing
+   end
 
-function formatter!(axis, low, high)
-    f = matplotlib.ticker.ScalarFormatter()
-    f.set_powerlimits((low, high))
-    axis.set_major_formatter(f)
-    return nothing
-end
+   """
+   $(TYPEDSIGNATURES)
+   """
+   function formatter!(axis, low, high)
+      f = matplotlib.ticker.ScalarFormatter()
+      f.set_powerlimits((low, high))
+      axis.set_major_formatter(f)
+      return nothing
+   end
+end #module
\ No newline at end of file
diff --git a/src/space.jl b/src/space.jl
index d6d48e2..c29dcca 100644
--- a/src/space.jl
+++ b/src/space.jl
@@ -1,79 +1,129 @@
 """
-numpy like linspace.
-
-$(SIGNATURES)
+$(TYPEDSIGNATURES)
+numpy like linspace.\\
+Args:\\
+- [Req] x0::Real    (start point)\\
+- [Req] xn::Real    (end point)\\
+- [Req] n::Integer  (# of points in stop-start range.)\\
+Note:\\
+^ endpoint false.\\
+^ useful for generating lattice of points in periodic box where the right hand edge is the box imaginary boundary.
 """
 function linspace(x0::Real, xn::Real, n::Integer)
-    return LinRange(x0, xn, n + 1)[1:end-1]
+   return LinRange(x0, xn, n + 1)[1:end-1]
 end
-
 """
-Mesh in cartesian coordinates.
+$(TYPEDSIGNATURES)
+"""
+function logspace(x0::Real, xn::Real, n::Integer;base=10)
+   return base.^(range(log(base,x0),stop=log(base,xn),length=n+1)[1:end-1])
+end
+#####################
 
-$(SIGNATURES)
+"""
+Mesh in cartesian coordinates by explictly passing x (row) and y (col) range vectors.
+Ouput is produced in classic (x,y) format with (id = row+nx*col) rule from left bottom corner
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] x::AbstractVector{T} (x coordinates range, row)\\
+[Req] y::AbstractVector{T} (y coordinates range, col. Use y=x if square box)\\
+Examples:\\
+cartesian_mesh(0:1:10,0:1:10)
+"""
+function cartesian_mesh(x::AbstractVector{T},y::AbstractVector{T}) where {T<:Real}
+   return (repeat(x, outer=length(y)), repeat(y, inner=length(x)))
+end
+"""
+Mesh in cartesian coordinates
+$(TYPEDSIGNATURES)
+Args:\\
+[Req] nx::Integer  (# of points in x axis)\\
+[Req] lx::Real     (length of the box in x axis)\\
+[Opt] ny = nx      (# of points in y axis)\\
+[Opt] ly = lx      (length of the box in y axis)\\
+[Opt] x0 = 0       (ofset coordinates in x axis)\\
+[Opt] y0 = 0       (ofset coordinates in y axis)\\
+Note:\\
+@ endpoint false. coordinates ranges x=[-lx/2, lx/2-1], y=[-ly/2, ly/2-1]\\
+@ useful for generating lattice of points in periodic box where the right hand edge is the box boundary.\\
+@ output is tuple follwoing (id = row+nx*col) rule from left bottom corner
 """
 function cartesian_mesh(nx::Integer, lx::Real; ny = nx, ly = lx, x0 = 0, y0 = 0)
-    x = linspace(-lx / 2, lx / 2, nx) .- x0
-    y = linspace(-ly / 2, ly / 2, ny)' .- y0
-    return x, y
+   x = linspace(-lx / 2, lx / 2, nx) .- x0
+   y = (linspace(-ly / 2, ly / 2, ny) .- y0)
+      return  x,y'
 end
+#####################
 
 """
-Mesh in polar coordinates.
-
-$(SIGNATURES)
+$(TYPEDSIGNATURES)
+Mesh in polar coordinates.\\
+Args:\\
+[Req] nr::Integer (# of points in radial direction)
+[Req] lr::Real;   (length of the radius)
+[Opt] nθ = nr     (# of points in θ direction)
+[Opt] lθ = 2π     (length of θ direction)
+[Opt] r0 = 0      (ofset in radial direction; not implemented yet)
+[Opt] θ0 = 0      (ofset in theta direction; not implemented yet)
 """
 function polar_mesh(nr::Integer, lr::Real; nθ = nr, lθ = 2π, r0 = 0, θ0 = 0)
-    r = linspace(r0, lr, nx)
-    θ = linspace(θ0, lθ, nθ)'
-    return r, θ
+   r = linspace(r0, lr, nr)
+   θ = linspace(θ0, lθ, nθ)
+      return  r, θ'
 end
 
 """
-Convert cartesian mesh to polar mesh.
-
-$(SIGNATURES)
+$(TYPEDSIGNATURES)
+Convert cartesian mesh to polar mesh.\\
+Args:\\
+[Req] x::AbstractVector{T} (explicit x coordinates in classic form, row)\\
+[Req] y::AbstractVector{T} (explicit y coordinates in classic form, col)\\
+[Opt] classic::Bool=false\\
+Examples:\\
+cartesian_mesh(0:1:10,0:1:10)\\
 """
-function polar_mesh(x::AbstractVector{T}, y::AbstractMatrix{T}) where {T<:Real}
-    r = hypot.(x, y)
-    θ = atan.(y, x)
-    return r, θ
+function polar_mesh(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Real}
+   r = hypot.(x, y)
+   θ = (atan.(y, x))
+   return  r, θ
 end
+#####################
 
-@inline function periodic_distance(x::T, l::T, hl::T) where {T<:Real}
-    @. ifelse(abs(x) > hl, l - x, x)
-end
-
-@inline function periodic_distance(x::T, l::T) where {T<:Real}
-    return periodic_distance(x, l, 0.5 * l)
+"""
+$(TYPEDSIGNATURES)
+Picks logspaced elements from given 1D array\\
+Args:\\
+[Req] x::AbstractVector;   (input vector of ordered data)\\
+[Opt] base = 1.2  (choose base )
+"""
+function logspacedPick(x::AbstractVector{T}; base = 1.2) where {T<:Real}
+   N = log(length(x)) / log(base)
+   n = unique(floor.(Int, base .^ (1:N)))
+   return x[n],n
 end
-
 """
-Convert periodic data to an aperiodic one.
-
-$(SIGNATURES)
-"""
-function aperiodic(x::Vector{T}, p = 2π) where {T<:Real}
-    y = zero(x)
-    x0, Δ, Δp = 0, 0, 0.9 * p
-    for (i, xi) in enumerate(x)
-        Δx = xi - x0
-        if (abs(Δx) >= Δp)
-            Δ = Δ - sign(Δx) * p
-        end
-        y[i] = xi + Δ
-        x0 = xi
-    end
-    return y
+$(TYPEDSIGNATURES)
+Picks logspaced elements from given 1D array but modifies the original input data.
+"""
+function logspacedPick!(x::AbstractVector{T}; base = 1.2) where {T<:Real}
+   sort!(x) #mutating orginal input vector
+   N = log(length(x)) / log(base)
+   n = unique(floor.(Int, base .^ (1:N)))
+   return x[n],n
 end
-
 """
-Pick logspaced elements from an array.
-
-$(SIGNATURES)
+$(TYPEDSIGNATURES)
+Picks logspaced radial elements from given 2D array after sorting the radius data.\\
+Output is in radius^2 format\\
+Args:\\
+[Req] x::AbstractVector;   (input vector x cords; explicit form)\\
+[Req] y::AbstractVector;   (input vector y cords; explicit form)\\
+[Opt] base = 1.2  (choose base)\\
 """
-function logspaced(x::AbstractVector; base = 1.2)
-    N = log(length(x)) / log(base)
-    n = unique(floor.(Int, base .^ (1:N)))
-    return x[n]
+function logspacedPick(x::AbstractVector{T},y::AbstractVector{T}; base = 1.2) where {T<:Real}
+   # r = @. (x^2 + y^2)^(0.5)
+   r = @. (x^2 + y^2)
+   sort!(r)
+   return logspacedPick(r,base=base)
 end
+#####################
\ No newline at end of file
diff --git a/src/spectrum.jl b/src/spectrum.jl
index 193c083..787dc1f 100644
--- a/src/spectrum.jl
+++ b/src/spectrum.jl
@@ -4,134 +4,181 @@ using SpecialFunctions: besselj0
 abstract type StructureFactor end
 abstract type Correlator end
 
+# define concrete data types of abstract data types.
 struct IsotropicStructureFactor <: StructureFactor
-    D  :: Integer
-    L  :: Float64
-    Δk :: Float64
-    k  :: Vector{Float64}
-    sk :: Vector{Float64}
+   D  :: Integer #dimensions
+   L  :: Float64
+   Δk :: Float64
+   k  :: Vector{Float64}
+   sk :: Vector{Float64}
 end
 
+# define constructors for the structs.
 function IsotropicStructureFactor(N, D; L = 2π)
-    Δk = 2π / L
-    k  = @. Δk * LinRange(0, N, N + 1)
-    sk = zeros(N + 1)
-    return IsotropicStructureFactor(D, L, Δk, k, sk)
+   Δk = 2π / L
+   k  = @. Δk * LinRange(0, N, N + 1)
+   sk = zeros(N + 1)
+   return IsotropicStructureFactor(D, L, Δk, k, sk)
 end
 
+# define concrete data types of abstract data types.
 struct IsotropicAutoCorrelator <: Correlator
-    D  :: Integer
-    r  :: Vector{Float64}
-    cr :: Vector{Float64}
+   D  :: Integer
+   r  :: Vector{Float64}
+   cr :: Vector{Float64}
 end
 
+# define constructors for the structs.
 function IsotropicAutoCorrelator(N, D; L = π)
-    r  = LinRange(0, L, N + 1)
-    cr = zeros(N + 1)
-    return IsotropicAutoCorrelator(D, r, cr)
+   r  = LinRange(0, L, N + 1)
+   cr = zeros(N + 1)
+   return IsotropicAutoCorrelator(D, r, cr)
 end
 
 # Returns the magnitude of the largest wavenumber for the range kiter.
 # (r)fftfreq omits the largest possible wavenumber so we add one.
 @inline function _kmax(kiter)
-    kmax = maximum(kiter) .+ 1
-    return round(Int, hypot(kmax...)) + 1
+   kmax = maximum(kiter) .+ 1
+   return round(Int, hypot(kmax...)) + 1
 end
 
+# function calculate kspace points in (1D,2D,3D) given if the input data is real or complex
 function _kgrid(nx::T, xfftfreq::Function) where {T<:Integer}
-    return Iterators.product(xfftfreq(nx))
+   return Iterators.product(xfftfreq(nx))
 end
 
 function _kgrid(nx::T, ny::T, xfftfreq::Function) where {T<:Integer}
-    return Iterators.product(xfftfreq(nx), fftfreq(ny, ny))
+   return Iterators.product(xfftfreq(nx), fftfreq(ny, ny))
 end
 
 function _kgrid(nx::T, ny::T, nz::T, xfftfreq::Function) where {T<:Integer}
-    return Iterators.product(xfftfreq(nx), fftfreq(ny, ny), fftfreq(nz, nz))
+   return Iterators.product(xfftfreq(nx), fftfreq(ny, ny), fftfreq(nz, nz))
 end
 
 @inline function _kbin(k::Tuple{Vararg{T}}) where {T<:Real}
-    return round(Int, hypot(k...)) + 1
+   return round(Int, hypot(k...)) + 1
 end
 
 # Computes the structure factor S(k) given the fourier transform uk = ⨏(u) of the data.
 # TODO: Write documentation.
 # The internals differ slightly for a real u and a complex u.
 
+# function to scale k₀ and kₙ if the data is real
 function _scale(uk::Array{ComplexF64}, factor)
-    uk0 = selectdim(uk, 1, 1)
-    ukN = selectdim(uk, 1, size(uk)[1])
-    @. uk0 = sqrt(0.5) * uk0
-    @. ukN = sqrt(0.5) * ukN
+   uk0 = selectdim(uk, 1, 1)
+   ukN = selectdim(uk, 1, size(uk)[1])
+   @. uk0 = factor * uk0
+   @. ukN = factor * ukN
 end
 
+# function to calculate structure_factor given FFT of data uk
+# calculating s(k) = |uₖ|² and binning in k values
 function structure_factor!(S::IsotropicStructureFactor, kiter::Iterators.ProductIterator, uk::Array{ComplexF64}, norm; isreal = true, preserve = true)
-    isreal ? _scale(uk, sqrt(0.5)) : nothing
+   isreal ? _scale(uk, sqrt(0.5)) : nothing
 
-    @inbounds for (i, k) in enumerate(kiter)
-        kb = _kbin(k)
-        S.sk[kb] = S.sk[kb] + abs2(uk[i])
-    end
-    @. S.sk = norm * S.sk
-    (isreal & preserve) ? _scale(uk, 1/sqrt(0.5)) : nothing
-    return nothing
+   @inbounds for (i, k) in enumerate(kiter)
+      kb = _kbin(k)
+      S.sk[kb] = S.sk[kb] + abs2(uk[i])
+   end
+   # @. S.sk = norm * S.sk
+   (isreal & preserve) ? _scale(uk, 1/sqrt(0.5)) : nothing
+   return nothing
 end
 
 function budget!( S::IsotropicStructureFactor, kiter::Iterators.ProductIterator, uk::Array{ComplexF64}, vk::Array{ComplexF64}, norm; isreal = true, preserve = true)
-    isreal ? _scale(uk, 0.5) : nothing
-    @inbounds for (i, k) in enumerate(kiter)
-        kb = _kbin(k)
-        S.sk[kb] = S.sk[kb] + real(conj(uk[i])*vk[i])
-    end
-    @. S.sk = norm * S.sk
-    (isreal & preserve) ? _scale(uk, 1/0.5) : nothing
-    return nothing
+   isreal ? _scale(uk, 0.5) : nothing
+   @inbounds for (i, k) in enumerate(kiter)
+      kb = _kbin(k)
+      S.sk[kb] = S.sk[kb] + real(conj(uk[i])*vk[i])
+   end
+   @. S.sk = norm * S.sk
+   (isreal & preserve) ? _scale(uk, 1/0.5) : nothing
+   return nothing
 end
 
+# function calculate kspace points depending upon if the input data is real or complex
 _xfftfreq(nx, ::Val{true})  = rfftfreq(2 * (nx - 1), 2 * (nx - 1))
 _xfftfreq(nx, ::Val{false}) = fftfreq(nx, nx)
 
+# function calculate normalization constant depending upon if the input data is real or complex
 _normalization(N, ::Val{true})  = 1.0 / (2 * (N[1] - 1) * prod(N[2:end]))^2
 _normalization(N, ::Val{false}) = 1.0 / (2.0 * prod(N)^2)
 
+# set up kgrid, normalization constant and initialize IsotropicStructureFactor (S)
 function _setup(N::Tuple{Vararg{Int}}; isreal = true, L = 2π)
-    kiter = _kgrid(N..., ((nx) -> _xfftfreq(nx, Val(isreal))))
-    S = IsotropicStructureFactor(_kmax(kiter), length(N); L = L)
-    return kiter, _normalization(N, Val(isreal)), S
-end
-
+   kiter = _kgrid(N..., ((nx) -> _xfftfreq(nx, Val(isreal))))
+   S = IsotropicStructureFactor(_kmax(kiter), length(N); L = L)
+   return kiter, _normalization(N, Val(isreal)), S
+end
+
+# main functions
+"""
+$(TYPEDSIGNATURES) 
+calculate isotropic structure factor for given `uk ≡ fft(data)` \\
+Args:\\
+
+   - `uk::Array{ComplexF64}` [Fourier transformed Array of input data. 1D, 2D or 3D arrays]\\
+   - `preserve = true` [will be depricated. if `true` prevents mutation of input data]\\
+   - `isreal = true` [read the note below]\\
+   - `kwargs...` [accepts, `L` size of the box, default is `2π`]\\
+
+!!! note
+    - use `true`  if `data` is real and `uk = rrft(data)`\\
+    - use `false` if `data` is real but `uk = fft(data)`, this use full range of fourier modes and expensive.\\
+    - use `false` if `data` is complex and `uk = fft(data)`\\
+"""
 function isotropic_structure_factor(uk::Array{ComplexF64}; preserve = true, isreal = true, kwargs...)
-    kiter, norm, S = _setup(size(uk); isreal, kwargs...)
-    structure_factor!(S, kiter, uk, norm; isreal, preserve)
-    return S
+   kiter, norm, S = _setup(size(uk); isreal, kwargs...)
+   structure_factor!(S, kiter, uk, norm; isreal, preserve)
+   @. S.sk = norm * S.sk
+   return S
 end
+"""
+$(TYPEDSIGNATURES) 
+calculate isotropic structure factor for given `uk = fft(data)` \\
+Args:\\
+
+   - `uk::Tuple{Vararg{Array{ComplexF64}}}`\\
 
+Takes tuple of Fourier transformed Array of input data. (1D, 2D or 3D arrays). Useful for a vectorial field.\\
+effectively Calculates `|uₖ|² + |vₖ|²` where `u` and `v` are components of vectorial field.
+"""
 function isotropic_structure_factor(uk::Tuple{Vararg{Array{ComplexF64}}}; preserve = true, isreal = true, kwargs...)
-    kiter, norm, S = _setup(size(uk[1]); isreal, kwargs...)
-    for uki in uk
-        structure_factor!(S, kiter, uki, norm; isreal, preserve)
-    end
-    return S
+   kiter, norm, S = _setup(size(uk[1]); isreal, kwargs...)
+   for uki in uk
+      structure_factor!(S, kiter, uki, norm; isreal, preserve)
+   end
+   @. S.sk = norm * S.sk
+   return S
 end
 
 function budget(uk::Array{ComplexF64}, vk::Array{ComplexF64}; preserve = true, isreal = true, kwargs...)
-    size(uk) == size(vk) || throw(DimensionMismatch(
-        "size(uk) = $(size(uk)) should be identical to size(vk) = $(size(vk))"
-    ))
-    kiter, norm, S = _setup(size(uk); isreal, kwargs...)
-    budget!(S, kiter, uk, vk, norm; isreal, preserve)
-    return S
+   size(uk) == size(vk) || throw(DimensionMismatch(
+      "size(uk) = $(size(uk)) should be identical to size(vk) = $(size(vk))"
+   ))
+   kiter, norm, S = _setup(size(uk); isreal, kwargs...)
+   budget!(S, kiter, uk, vk, norm; isreal, preserve)
+   return S
 end
 
 @inline _expintegral(kr, ::Val{1}) = cos(kr)
 @inline _expintegral(kr, ::Val{2}) = besselj0(kr)
 @inline _expintegral(kr, ::Val{3}) = sinc(kr/π)
 
+"""
+$(TYPEDSIGNATURES) 
+Calculates isotropic correlation function `C(r) = ⟨u(0)u(r)⟩` given isotropic structure function `S(k)`\\
+Args:\\
+
+   - `S::IsotropicStructureFactor` [StructureFactor object from [`isotropic_structure_factor`](@ref) function]\\
+   - `N = 128` [number of discrete `r` points for `C(r)`, keep it below the box size, Nyquist theorem]
+
+"""
 function correlation(S::IsotropicStructureFactor; N = 128)
-    C = IsotropicAutoCorrelator(N, S.D; L = S.L / 2)
-    @inbounds for (i, r) in enumerate(C.r)
-        C.cr[i] = sum((@. S.sk * _expintegral(S.k * r, Val(S.D))))
-    end
-    @. C.cr = C.cr / C.cr[1]
-    return C
+   C = IsotropicAutoCorrelator(N, S.D; L = S.L / 2)
+   @inbounds for (i, r) in enumerate(C.r)
+      C.cr[i] = sum((@. S.sk * _expintegral(S.k * r, Val(S.D))))
+   end
+   @. C.cr = C.cr / C.cr[1]
+   return C
 end