diff --git a/src/pair-correlation.jl b/src/pair-correlation.jl
index 56f4e92..d92a26d 100644
--- a/src/pair-correlation.jl
+++ b/src/pair-correlation.jl
@@ -89,36 +89,37 @@ struct DiscretePairCorrelation{Dim} <: PairCorrelation
         end
         
         # Predict the number density from the restriction that is due to the link with a probability distribution, and assuming homogeneous distribution of particles. 
+        if !isempty(r)
+            σs = trapezoidal_scheme(r)
 
-        σs = trapezoidal_scheme(r)
+            sumg = sum(g .* r .^ (Dim - 1) .* σs)
 
-        sumg = sum(g .* r .^ (Dim - 1) .* σs)
+            number_density_predict = 1 / (2*(Dim - 1) * π * (r[end] ^ Dim / Dim - sumg))
 
-        number_density_predict = 1 / (2*(Dim - 1) * π * (r[end] ^ Dim / Dim - sumg))
+            if number_density == -1.0
+                number_density = number_density_predict
+            end
 
-        if number_density == -1.0
-            number_density = number_density_predict
-        end
-
-        number_density_error = abs(number_density / number_density_predict - 1)
+            number_density_error = abs(number_density / number_density_predict - 1)
         
-        # This formula seems to sensitive to numerical errors.
-        # if number_density_error > tol
-        #     println("The number density that was specified was $(number_density), whereas the calculated number density was $(number_density_predict), a reltive error of $(number_density_error)")
-        # end
+            # This formula seems to sensitive to numerical errors.
+            # if number_density_error > tol
+            #     println("The number density that was specified was $(number_density), whereas the calculated number density was $(number_density_predict), a reltive error of $(number_density_error)")
+            # end
         
-        if rescale
-            println("rescaling this pair correlation according to the number density provided.")
-            # this rescaling is due to a restriction imposed by the connection to probability distributions
+            if rescale
+                println("rescaling this pair correlation according to the number density provided.")
+                # this rescaling is due to a restriction imposed by the connection to probability distributions
 
-            b = r[end] ^ Dim / Dim - 1 / (2*(Dim - 1) * π * number_density)
-            α = b / sumg
+                b = r[end] ^ Dim / Dim - 1 / (2*(Dim - 1) * π * number_density)
+                α = b / sumg
 
-            if α < 0
-                @error "rescaling has failed and has given a negative scale for the pair-correlation"
-            end     
+                if α < 0
+                    @error "rescaling has failed and has given a negative scale for the pair-correlation"
+                end     
 
-            g = α .* g
+                g = α .* g
+            end
         end
 
         new{Dim}(r,g,minimal_distance,number_density)
diff --git a/test/particulate-from-structure.jl b/test/particulate-from-structure.jl
index cb861ca..1fc783d 100644
--- a/test/particulate-from-structure.jl
+++ b/test/particulate-from-structure.jl
@@ -3,12 +3,77 @@
 using ParticleCorrelations
 using Test, Statistics
 
-using Optim
 using Plots
 
+# let us first aim to recover a MonteCarlo pair-correlation from a crystalline arrangement
+
+# choose the spatial dimension
+dim = 2
+
+# choose the medium for the particles. Currently only one type
+medium = HardMedium{dim}()
+
+# choose the shapes of the particles
+
+pairtype = MonteCarloPairCorrelation(dim; 
+    rtol = 1e-3, 
+    maxlength = 100, 
+    iterations = 10, 
+    numberofparticles = 3000
+)
+
+# choose the medium for the particles. Currently only one type
+medium = HardMedium{dim}()
+
+# Choose the species, which represents a collection of one type of particle
+radius = 0.5
+
+s = Specie(
+    medium,
+    Sphere(dim, radius),
+    volume_fraction = 0.15,
+    separation_ratio = 1.0 # minimal distance from this particle = r * (separation_ratio - 1.0) 
+);
+
+rs = 0.2:0.4:8.0
+pair = pair_correlation(s, pairtype, rs)
+
+
+# If you have the Plots package
+plot(pair.r, pair.g)
+
+ks = 1.0:0.5:20.0
+sfactor = structure_factor(pair, ks)
+
+plot(sfactor.k, sfactor.S)
+
+
+# Define two regions R1 and R2 to place the particles
+reg1 = Box([[-40.0, -40.0],[40.0, 40.0]]);
+reg2 = Box([[-42.0, -42.0],[42.0, 42.0]]);
+
+Dim = 2 
+
+# box = reg2
+function periodic_particles(box::Box{T,Dim}, specie::Specie{Dim}) where {T, Dim}
+    
+    Id = Matrix(I, Dim, Dim)
+    vs = [Id[:,i] for i = 1:Dim]
+
+    n = number_density(specie)
+
+    cell_length = (1 / n)^(1/Dim) / 2
+
+    d1 = origin(box) - box.dimensions ./ (2)
+    d2 = origin(box) + box.dimensions ./ (2)
+    
+end    
+
+using Optim
 
-a = 1.0
-rs = 0.0:0.01:8.0
+function g!(G, S)
+    G[1] = -2.0 * (1.0 - S[1]) 
+end
 
 g = 1.0 .+ exp.(-(rs .- 4a).^2)