diff --git a/Project.toml b/Project.toml
index d97de214..48f153d1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,7 +8,9 @@ Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 PlotThemes = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
@@ -26,8 +28,8 @@ DocumenterTools = "0.1"
 Interpolations = "0.13"
 NumericalIntegration = "0.2"
 Parameters = "0.12"
-PlotlyBase = "0.4"
 PlotThemes = "2.0"
+PlotlyBase = "0.4"
 PyCall = "1.92"
 QuadGK = "2.4"
 SciPy = "0.1"
diff --git a/src/AngularCoefficients.jl b/src/AngularCoefficients.jl
index 90ca0405..53f363d6 100644
--- a/src/AngularCoefficients.jl
+++ b/src/AngularCoefficients.jl
@@ -1,18 +1,24 @@
+abstract type IntegrationMethod end
+
+struct NumericalIntegrationSimpson <: IntegrationMethod end
+struct CustomTrapz <: IntegrationMethod end
+
+
 """
     ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients,
     WeightFunctionA::GCWeightFunction, WeightFunctionB::GCWeightFunction,
     BackgroundQuantities::BackgroundQuantities,
     w0waCDMCosmology::w0waCDMCosmology, CosmologicalGrid::CosmologicalGrid,
-    PowerSpectrum::PowerSpectrum)
+    PowerSpectrum::PowerSpectrum, ::NumericalIntegrationSimpson)
 
 This function evaluates the Angular Coefficients for all tomographic bins and
 multipole values.
 """
 function  ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients,
-    WeightFunctionA::GCWeightFunction, WeightFunctionB::GCWeightFunction,
+    WeightFunctionA::WeightFunction, WeightFunctionB::WeightFunction,
     BackgroundQuantities::BackgroundQuantities,
     w0waCDMCosmology::w0waCDMCosmology, CosmologicalGrid::CosmologicalGrid,
-    PowerSpectrum::PowerSpectrum)
+    PowerSpectrum::PowerSpectrum, ::NumericalIntegrationSimpson)
     c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                        #physical constants involved in calculations
     for idx_a in 1:size(AngularCoefficients.AngularCoefficientsArray, 2)
@@ -27,7 +33,36 @@ function  ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients,
                 AngularCoefficients.AngularCoefficientsArray[idx_l, idx_a,
                 idx_b] = NumericalIntegration.integrate(CosmologicalGrid.ZArray,
                 Integrand, SimpsonEvenFast())
+                AngularCoefficients.AngularCoefficientsArray[idx_l, idx_b,
+                idx_a] = AngularCoefficients.AngularCoefficientsArray[idx_l,
+                idx_a, idx_b]
             end
         end
     end
 end
+
+function  ComputeAngularCoefficients(AngularCoefficients::AngularCoefficients,
+    WeightFunctionA::WeightFunction, WeightFunctionB::WeightFunction,
+    BackgroundQuantities::BackgroundQuantities,
+    w0waCDMCosmology::w0waCDMCosmology, CosmologicalGrid::CosmologicalGrid,
+    PowerSpectrum::PowerSpectrum, ::CustomTrapz)
+    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
+                       #physical constants involved in calculations
+    Integrand = zeros(size(AngularCoefficients.AngularCoefficientsArray,1),
+    size(AngularCoefficients.AngularCoefficientsArray,2),
+    size(AngularCoefficients.AngularCoefficientsArray,3))
+    @avx for i ∈ axes(AngularCoefficients.AngularCoefficientsArray,2),
+        j ∈ axes(AngularCoefficients.AngularCoefficientsArray,3),
+        l ∈ axes(AngularCoefficients.AngularCoefficientsArray,1)
+        for z ∈ axes(CosmologicalGrid.ZArray,1)
+            Integrand[l,i,j] += c_0 *
+            WeightFunctionA.WeightFunctionArray[i, z] *
+            WeightFunctionB.WeightFunctionArray[j, z] /
+            (BackgroundQuantities.HZArray[z] *
+            BackgroundQuantities.rZArray[z]^2) *
+            PowerSpectrum.InterpolatedPowerSpectrum[l,z]
+        end
+    end
+    Integrand .*= (last(CosmologicalGrid.ZArray)-first(CosmologicalGrid.ZArray))/(length(CosmologicalGrid.ZArray)-1)
+    AngularCoefficients.AngularCoefficientsArray = Integrand
+end
diff --git a/src/CosmoCentral.jl b/src/CosmoCentral.jl
index dfdd5720..8c8b8bee 100644
--- a/src/CosmoCentral.jl
+++ b/src/CosmoCentral.jl
@@ -4,8 +4,11 @@ using QuadGK
 using Base: @kwdef
 using Parameters
 using NumericalIntegration
+using Interpolations
 using PyCall
 using Dierckx
+using HDF5
+using LoopVectorization
 numpy = pyimport("numpy")
 classy = pyimport("classy")
 
@@ -18,5 +21,6 @@ include("WeightFunctions.jl")
 include("BoltzmannSolver.jl")
 include("PowerSpectrum.jl")
 include("AngularCoefficients.jl")
+include("FileIO.jl")
 
 end # module
diff --git a/src/CosmologicalStructures.jl b/src/CosmologicalStructures.jl
index 567e276b..85531604 100644
--- a/src/CosmologicalStructures.jl
+++ b/src/CosmologicalStructures.jl
@@ -13,6 +13,7 @@ abstract type ConvolvedDensity <: AsbtractDensity end
 abstract type InstrumentResponse end
 abstract type WeightFunction end
 abstract type GCWeightFunction <: WeightFunction end
+abstract type WLWeightFunction <: WeightFunction end
 abstract type PowerSpectrum end
 abstract type AngularCoefficients end
 
@@ -143,7 +144,7 @@ The parameters contained in this struct are
 @kwdef mutable struct AnalitycalDensityStruct <: AnalitycalDensity
     Z0::Float64 = 0.9/sqrt(2.)
     ZMin::Float64 = 0.001
-    ZMax::Float64 = 2.5
+    ZMax::Float64 = 4.0
     SurfaceDensity::Float64 = 30.
     Normalization::Float64 = 1.
 end
@@ -211,6 +212,18 @@ for all tomographic bins and redshift values in the [`CosmologicalGridStruct`](@
     WeightFunctionArray::AbstractArray{Float64, 2}
 end
 
+"""
+    WLWeightFunctionStruct()
+
+This struct contains the array with the Weak Lensing Weight function values
+for all tomographic bins and redshift values in the [`CosmologicalGridStruct`](@ref).
+"""
+@kwdef mutable struct WLWeightFunctionStruct <: WLWeightFunction
+    WeightFunctionArray::AbstractArray{Float64, 2}
+    LensingEfficiencyArray::AbstractArray{Float64, 2}
+end
+
+
 """
     PowerSpectrumStruct()
 
@@ -235,3 +248,9 @@ This struct contains the array with the Angular Coefficients.
 @kwdef mutable struct AngularCoefficientsStruct <: AngularCoefficients
     AngularCoefficientsArray::AbstractArray{Float64, 3} = zeros(2991, 10, 10)
 end
+
+abstract type InterpolationMethod end
+
+struct RectBivSplineDierckx <: InterpolationMethod end
+struct GriddedLinear <: InterpolationMethod end
+struct BSplineCubic <: InterpolationMethod end
diff --git a/src/FileIO.jl b/src/FileIO.jl
new file mode 100644
index 00000000..a4981153
--- /dev/null
+++ b/src/FileIO.jl
@@ -0,0 +1,90 @@
+function WriteAngularCoefficients(AngularCoefficients::AngularCoefficients,
+    CosmologicalGrid::CosmologicalGrid, WeightFunction::WeightFunction,
+    Bias::Bias, ConvolvedDensity::ConvolvedDensity, Filename::String)
+    h5write(Filename*".h5", "cls/PhotometricGalaxy_PhotometricGalaxy/c_lij",
+    AngularCoefficients.AngularCoefficientsArray)
+    h5write(Filename*".h5",
+    "weight_functions/PhotometricGalaxy/bias/b_i_z",
+    Bias.BiasArray)
+    h5write(Filename*".h5",
+    "weight_functions/PhotometricGalaxy/density/norm_density_iz",
+    ConvolvedDensity.DensityNormalizationArray)
+end
+
+function ReadAngularCoefficients(Filename::String)
+    Filename *= ".h5"
+    file = HDF5.h5open(Filename, "r")
+    c_lij =
+    HDF5.read(file["cls"]["PhotometricGalaxy_PhotometricGalaxy"]["c_lij"])
+    AngularCoefficients = AngularCoefficientsStruct(AngularCoefficientsArray =
+    c_lij)
+    return AngularCoefficients
+end
+
+function WritePowerSpectrumBackground(PowerSpectrum::PowerSpectrum,
+    BackgroundQuantities::BackgroundQuantities,
+    CosmologicalGrid::CosmologicalGrid, Filename::String)
+    h5write(Filename*".h5",
+    "cosmology/comoving_distance_array",
+    BackgroundQuantities.rZArray)
+    h5write(Filename*".h5",
+    "cosmology/hubble_array",
+    BackgroundQuantities.HZArray)
+    h5write(Filename*".h5",
+    "cosmology/z_grid",
+    CosmologicalGrid.ZArray)
+    h5write(Filename*".h5",
+    "power_spectrum/z_grid",
+    CosmologicalGrid.ZArray)
+    h5write(Filename*".h5",
+    "power_spectrum/k_grid",
+    CosmologicalGrid.KArray)
+    h5write(Filename*".h5",
+    "power_spectrum/lin_p_mm_k_z",
+    PowerSpectrum.PowerSpectrumLinArray)
+    h5write(Filename*".h5",
+    "power_spectrum/nonlin_p_mm_k_z",
+    PowerSpectrum.PowerSpectrumNonlinArray)
+end
+
+function ReadPowerSpectrumBackground(Filename::String,
+    MultipolesArray::Vector{Float64})
+    Filename *= ".h5"
+    file = HDF5.h5open(Filename, "r")
+    nonlin_p_mm_k_z = HDF5.read(file["power_spectrum"]["nonlin_p_mm_k_z"])
+    lin_p_mm_k_z    = HDF5.read(file["power_spectrum"]["lin_p_mm_k_z"])
+    z_grid = HDF5.read(file["power_spectrum"]["z_grid"])
+    k_grid = HDF5.read(file["power_spectrum"]["k_grid"])
+    comoving_distance_array = HDF5.read(file["cosmology"]["comoving_distance_array"])
+    hubble_array = HDF5.read(file["cosmology"]["hubble_array"])
+    BackgroundQuantities = BackgroundQuantitiesStruct(HZArray = hubble_array,
+    rZArray = comoving_distance_array)
+    CosmologicalGrid = CosmologicalGridStruct(ZArray = z_grid, KArray = k_grid,
+    MultipolesArray = MultipolesArray)
+    PowerSpectrum = PowerSpectrumStruct(PowerSpectrumLinArray = lin_p_mm_k_z,
+    PowerSpectrumNonlinArray = nonlin_p_mm_k_z,
+    InterpolatedPowerSpectrum = zeros(length(CosmologicalGrid.MultipolesArray), length(CosmologicalGrid.ZArray)))
+    return PowerSpectrum, BackgroundQuantities, CosmologicalGrid
+end
+
+function ReadPowerSpectrumBackgroundSeyfert(Filename::String,
+    MultipolesArray::Vector{Float64})
+    c_0 = 2.99792458e5
+    Filename *= ".h5"
+    file = HDF5.h5open(Filename, "r")
+    nonlin_p_mm_k_z = HDF5.read(file["power_spectrum"]["nonlin_p_mm_z_k"])
+    lin_p_mm_k_z    = HDF5.read(file["power_spectrum"]["lin_p_mm_z_k"])
+    z_grid = HDF5.read(file["power_spectrum"]["z_grid"])
+    k_grid = HDF5.read(file["power_spectrum"]["k_grid"])
+    dimensionless_comoving_distance_array = HDF5.read(file["cosmology"]["dimensionless_comoving_distance_array"])
+    dimensionless_hubble_array = HDF5.read(file["cosmology"]["dimensionless_hubble_array"])
+    BackgroundQuantities = BackgroundQuantitiesStruct(HZArray =
+    dimensionless_hubble_array*67,
+    rZArray = dimensionless_comoving_distance_array*c_0/67)
+    CosmologicalGrid = CosmologicalGridStruct(ZArray = z_grid, KArray = k_grid,
+    MultipolesArray = MultipolesArray)
+    PowerSpectrum = PowerSpectrumStruct(PowerSpectrumLinArray = lin_p_mm_k_z,
+    PowerSpectrumNonlinArray = nonlin_p_mm_k_z,
+    InterpolatedPowerSpectrum = zeros(length(CosmologicalGrid.MultipolesArray), length(CosmologicalGrid.ZArray)))
+    return PowerSpectrum, BackgroundQuantities, CosmologicalGrid
+end
diff --git a/src/PowerSpectrum.jl b/src/PowerSpectrum.jl
index 1257282f..afa5297b 100644
--- a/src/PowerSpectrum.jl
+++ b/src/PowerSpectrum.jl
@@ -24,7 +24,8 @@ This function interpolates the Power Spectrum on the ``k-z`` grid and evaluates
 it on the Limber grid.
 """
 function InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid::CosmologicalGrid,
-    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum)
+    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
+    ::RectBivSplineDierckx)
     InterpPmm = Dierckx.Spline2D(log10.(CosmologicalGrid.KArray),
     CosmologicalGrid.ZArray, log10.(PowerSpectrum.PowerSpectrumNonlinArray);
     kx=5, ky=5, s=0.0)
@@ -34,3 +35,36 @@ function InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid::CosmologicalGrid,
         CosmologicalGrid.ZArray))
     end
 end
+
+function InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
+    ::GriddedLinear)
+    InterpPmm = Interpolations.interpolate((log10.(CosmologicalGrid.KArray),
+    CosmologicalGrid.ZArray), log10.(PowerSpectrum.PowerSpectrumNonlinArray),
+    Gridded(Linear()))
+    InterpPmm = Interpolations.extrapolate(InterpPmm, Line())
+    for (idx_l, myl) in enumerate(CosmologicalGrid.MultipolesArray)
+        PowerSpectrum.InterpolatedPowerSpectrum[idx_l, :] =
+        10 .^(InterpPmm.(log10.(CosmologicalGrid.KLimberArray[idx_l, :]),
+        CosmologicalGrid.ZArray))
+    end
+end
+
+function InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities, PowerSpectrum::PowerSpectrum,
+    ::BSplineCubic)
+    x = LinRange(log10(first(CosmologicalGrid.KArray)),
+    log10(last(CosmologicalGrid.KArray)), length(CosmologicalGrid.KArray))
+    y = LinRange(first(CosmologicalGrid.ZArray), last(CosmologicalGrid.ZArray),
+    length(CosmologicalGrid.ZArray))
+    InterpPmm = Interpolations.interpolate(
+    log10.(PowerSpectrum.PowerSpectrumNonlinArray),
+    BSpline(Cubic(Line(OnGrid()))))
+    InterpPmm = scale(InterpPmm, x, y)
+    InterpPmm = Interpolations.extrapolate(InterpPmm, Line())
+    for (idx_l, myl) in enumerate(CosmologicalGrid.MultipolesArray)
+        PowerSpectrum.InterpolatedPowerSpectrum[idx_l, :] =
+        10 .^(InterpPmm.(log10.(CosmologicalGrid.KLimberArray[idx_l, :]),
+        CosmologicalGrid.ZArray))
+    end
+end
diff --git a/src/WeightFunctions.jl b/src/WeightFunctions.jl
index 789ca58f..094c1466 100644
--- a/src/WeightFunctions.jl
+++ b/src/WeightFunctions.jl
@@ -47,3 +47,83 @@ function ComputeWeightFunctionOverGrid(
         end
     end
 end
+
+"""
+    ComputeWeightFunction(z::Float64, i::Int64,
+        ConvolvedDensity::ConvolvedDensity,
+        w0waCDMCosmology::w0waCDMCosmology)
+
+This function returns the Galaxy Clustering Weight function for a given redshift
+``z`` and tomographic bin ``i``.
+"""
+function ComputeLensingEfficiency(z::Float64, i::Int64,
+    ConvolvedDensity::ConvolvedDensity, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, w0waCDMCosmology::w0waCDMCosmology,
+    CosmologicalGrid::CosmologicalGrid, WLWeightFunction::WLWeightFunction)
+    int, err = QuadGK.quadgk(x -> CosmoCentral.ComputeConvolvedDensityFunction(x, i,
+    ConvolvedDensity, AnalitycalDensity, InstrumentResponse)*
+    (1. - CosmoCentral.ComputeComovingDistance(z, w0waCDMCosmology)/
+    CosmoCentral.ComputeComovingDistance(x, w0waCDMCosmology)), z,
+    last(CosmologicalGrid.ZArray) , rtol=1e-12)
+    return int
+end
+
+"""
+    ComputeWeightFunctionOverGrid(
+    GCWeightFunction::GCWeightFunction, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, ConvolvedDensity::ConvolvedDensity,
+    Bias::Bias, CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities,
+    w0waCDMCosmology::w0waCDMCosmology)
+
+This function evaluates the Galaxy Clustering Weight Function over the ``z``
+grid and for all tomographic bins ``i``.
+"""
+function ComputeLensingEfficiencyOverGrid(
+    WLWeightFunction::WLWeightFunction, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, ConvolvedDensity::ConvolvedDensity,
+    CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities,
+    w0waCDMCosmology::w0waCDMCosmology)
+    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
+        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
+            WLWeightFunction.LensingEfficiencyArray[idx_ZBinArray, idx_ZArray] =
+            ComputeLensingEfficiency(CosmologicalGrid.ZArray[idx_ZArray],
+            idx_ZBinArray, ConvolvedDensity, AnalitycalDensity,
+            InstrumentResponse, w0waCDMCosmology, CosmologicalGrid,
+            WLWeightFunction)
+        end
+    end
+end
+
+function ComputeWeightFunction(z::Float64, i::Int64,
+    ConvolvedDensity::ConvolvedDensity, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, w0waCDMCosmology::w0waCDMCosmology,
+    CosmologicalGrid::CosmologicalGrid, WLWeightFunction::WLWeightFunction)
+    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
+                       #physical constants involved in calculations
+    return 1.5 * ComputeLensingEfficiency(z, i, ConvolvedDensity,
+    AnalitycalDensity,InstrumentResponse, w0waCDMCosmology, CosmologicalGrid,
+    WLWeightFunction) * (w0waCDMCosmology.H0/c_0)^2 * w0waCDMCosmology.ΩM *
+    (1. + z) * ComputeComovingDistance(z, w0waCDMCosmology)
+end
+
+
+function ComputeWeightFunctionOverGrid(
+    WLWeightFunction::WLWeightFunction, AnalitycalDensity::AnalitycalDensity,
+    InstrumentResponse::InstrumentResponse, ConvolvedDensity::ConvolvedDensity,
+    CosmologicalGrid::CosmologicalGrid,
+    BackgroundQuantities::BackgroundQuantities,
+    w0waCDMCosmology::w0waCDMCosmology)
+    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
+                       #physical constants involved in calculations
+    for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
+        for idx_ZArray in 1:length(CosmologicalGrid.ZArray)
+            WLWeightFunction.WeightFunctionArray[idx_ZBinArray, idx_ZArray] =
+            1.5 * (w0waCDMCosmology.H0/c_0)^2 * w0waCDMCosmology.ΩM *
+            (1. + CosmologicalGrid.ZArray[idx_ZArray]) *
+            BackgroundQuantities.rZArray[idx_ZArray] *
+            WLWeightFunction.LensingEfficiencyArray[idx_ZBinArray, idx_ZArray]
+        end
+    end
+end
diff --git a/test/Project.toml b/test/Project.toml
index b718e96c..ed5ff6a0 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,6 +1,9 @@
 [deps]
 Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 NumericalIntegration = "e7bfaba1-d571-5449-8927-abc22e82249b"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
diff --git a/test/cl.h5 b/test/cl.h5
new file mode 100644
index 00000000..d5e80010
Binary files /dev/null and b/test/cl.h5 differ
diff --git a/test/p_mm.h5 b/test/p_mm.h5
new file mode 100644
index 00000000..fc1cc0f7
Binary files /dev/null and b/test/p_mm.h5 differ
diff --git a/test/runtests.jl b/test/runtests.jl
index 7cc4e860..5cf82b40 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,16 +4,16 @@ using Test
 using QuadGK
 using NumericalIntegration
 using PyCall
+using Interpolations
 
 numpy = pyimport("numpy")
-classy = pyimport("classy")
 
 w0waCDMCosmology = CosmoCentral.w0waCDMCosmologyStruct()
 AnalitycalDensity = CosmoCentral.AnalitycalDensityStruct()
 InstrumentResponse = CosmoCentral.InstrumentResponseStruct()
 ConvolvedDensity = CosmoCentral.ConvolvedDensityStruct()
 CosmologicalGrid  = CosmoCentral.CosmologicalGridStruct(
-ZArray=Array(LinRange(0.001, 2.5, 300)))
+ZArray=Array(LinRange(0.001, 4.0, 500)))
 BackgroundQuantities = CosmoCentral.BackgroundQuantitiesStruct(HZArray=
 zeros(length(CosmologicalGrid.ZArray)),
 rZArray=zeros(length(CosmologicalGrid.ZArray)))
@@ -23,6 +23,12 @@ length(CosmologicalGrid.ZArray)))
 GCWeightFunction = CosmoCentral.GCWeightFunctionStruct(WeightFunctionArray =
 zeros(length(ConvolvedDensity.DensityNormalizationArray),
 length(CosmologicalGrid.ZArray)))
+WLWeightFunction = CosmoCentral.WLWeightFunctionStruct(WeightFunctionArray =
+zeros(length(ConvolvedDensity.DensityNormalizationArray),
+length(CosmologicalGrid.ZArray)),
+LensingEfficiencyArray = zeros(length(
+ConvolvedDensity.DensityNormalizationArray),
+length(CosmologicalGrid.ZArray)))
 AngularCoefficients = CosmoCentral.AngularCoefficientsStruct(
 AngularCoefficientsArray = zeros(length(CosmologicalGrid.MultipolesArray),
 length(GCWeightFunction.WeightFunctionArray[:, 1]),
@@ -119,47 +125,120 @@ end
 end
 
 @testset "Check the Weight function evaluation" begin
-    test_array = zeros(length(ConvolvedDensity.ZBinArray)-1,
+    test_gc = zeros(length(ConvolvedDensity.ZBinArray)-1,
+    length(CosmologicalGrid.ZArray))
+    test_le = zeros(length(ConvolvedDensity.ZBinArray)-1,
+    length(CosmologicalGrid.ZArray))
+    test_wl = zeros(length(ConvolvedDensity.ZBinArray)-1,
     length(CosmologicalGrid.ZArray))
     for (idxz, zvalue) in enumerate(CosmologicalGrid.ZArray)
         for idx_ZBinArray in 1:length(ConvolvedDensity.ZBinArray)-1
-            test_array[idx_ZBinArray, idxz] =
+            test_gc[idx_ZBinArray, idxz] =
+            CosmoCentral.ComputeWeightFunction(zvalue, idx_ZBinArray,
+            ConvolvedDensity, AnalitycalDensity,
+            InstrumentResponse, w0waCDMCosmology,
+            PiecewiseBias, GCWeightFunction)
+            test_le[idx_ZBinArray, idxz] =
+            CosmoCentral.ComputeLensingEfficiency(zvalue, idx_ZBinArray,
+            ConvolvedDensity, AnalitycalDensity, InstrumentResponse,
+            w0waCDMCosmology, CosmologicalGrid, WLWeightFunction)
+            test_wl[idx_ZBinArray, idxz] =
             CosmoCentral.ComputeWeightFunction(zvalue, idx_ZBinArray,
-        ConvolvedDensity, AnalitycalDensity,
-        InstrumentResponse, w0waCDMCosmology,
-        PiecewiseBias, GCWeightFunction)
+                ConvolvedDensity, AnalitycalDensity,
+                InstrumentResponse, w0waCDMCosmology, CosmologicalGrid,
+                WLWeightFunction)
         end
     end
-CosmoCentral.ComputeWeightFunctionOverGrid(GCWeightFunction, AnalitycalDensity, InstrumentResponse, ConvolvedDensity, PiecewiseBias, CosmologicalGrid, BackgroundQuantities, w0waCDMCosmology)
-@test isapprox(test_array, GCWeightFunction.WeightFunctionArray, atol=1e-9)
+    CosmoCentral.ComputeWeightFunctionOverGrid(GCWeightFunction,
+    AnalitycalDensity, InstrumentResponse, ConvolvedDensity, PiecewiseBias,
+    CosmologicalGrid, BackgroundQuantities, w0waCDMCosmology)
+    CosmoCentral.ComputeLensingEfficiencyOverGrid(
+    WLWeightFunction, AnalitycalDensity,
+    InstrumentResponse, ConvolvedDensity,
+    CosmologicalGrid,
+    BackgroundQuantities,
+    w0waCDMCosmology)
+    CosmoCentral.ComputeWeightFunctionOverGrid(WLWeightFunction,
+    AnalitycalDensity, InstrumentResponse, ConvolvedDensity, CosmologicalGrid,
+    BackgroundQuantities, w0waCDMCosmology)
+    @test isapprox(test_gc, GCWeightFunction.WeightFunctionArray, rtol=1e-9)
+    @test isapprox(test_wl, WLWeightFunction.WeightFunctionArray, rtol=1e-9)
+    @test isapprox(test_le, WLWeightFunction.LensingEfficiencyArray, rtol=1e-9)
 end
 
 
 @testset "Check the Power Spectrum evaluated over the Limber Grid" begin
+    MultipolesArray = Array(LinRange(10.5, 2999.5, 2990))
+    PowerSpectrum, BackgroundQuantitiesLoaded, CosmologicalGrid =
+    CosmoCentral.ReadPowerSpectrumBackground(
+    "/home/runner/work/CosmoCentral.jl/CosmoCentral.jl/test/p_mm",
+    MultipolesArray)
+    CosmoCentral.InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid,
+    BackgroundQuantitiesLoaded, PowerSpectrum, CosmoCentral.BSplineCubic())
     classyParams = CosmoCentral.Initializeclassy(w0waCDMCosmology)
-    PowerSpectrum = CosmoCentral.PowerSpectrumStruct(PowerSpectrumLinArray =
-    zeros(length(CosmologicalGrid.KArray), length(CosmologicalGrid.ZArray)),
-    PowerSpectrumNonlinArray = zeros(length(CosmologicalGrid.KArray),
-    length(CosmologicalGrid.ZArray)),
-    InterpolatedPowerSpectrum = zeros(length(CosmologicalGrid.MultipolesArray),
-    length(CosmologicalGrid.ZArray)))
-    CosmoCentral.ComputeLimberArray(CosmologicalGrid, BackgroundQuantities)
+    CosmoCentral.ComputeLimberArray(CosmologicalGrid,
+    BackgroundQuantitiesLoaded)
     test_k_limber = (CosmologicalGrid.MultipolesArray[1]+0.5) /
-    BackgroundQuantities.rZArray[1]
+    BackgroundQuantitiesLoaded.rZArray[1]
     @test test_k_limber == CosmologicalGrid.KLimberArray[1, 1]
     test_Omega_cdm = w0waCDMCosmology.ΩM-w0waCDMCosmology.ΩB-
     w0waCDMCosmology.Mν/(93.14*(w0waCDMCosmology.H0/100)^2)
     @test test_Omega_cdm == classyParams.classyParamsDict["Omega_cdm"]
-    #CosmoCentral.EvaluatePowerSpectrum(classyParams, CosmologicalGrid,
-    #PowerSpectrum)
-    #CosmoCentral.InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid,
-    #BackgroundQuantities, PowerSpectrum)
-    #cosmo = classy.Class()
-    #cosmo.set(classyParams.classyParamsDict)
-    #cosmo.compute()
-    #println(cosmo.pk(test_k_limber, CosmologicalGrid.ZArray[1]))
-    #test_power_spectrum = cosmo.pk(test_k_limber, CosmologicalGrid.ZArray[1])
-    #println(test_power_spectrum)
-    #@test isapprox(test_power_spectrum,
-    #PowerSpectrum.InterpolatedPowerSpectrum[1, 1], atol=1e-7)
+    x = LinRange(log10(first(CosmologicalGrid.KArray)),
+    log10(last(CosmologicalGrid.KArray)), length(CosmologicalGrid.KArray))
+    y = LinRange(first(CosmologicalGrid.ZArray), last(CosmologicalGrid.ZArray),
+    length(CosmologicalGrid.ZArray))
+    InterpPmm = Interpolations.interpolate(
+    log10.(PowerSpectrum.PowerSpectrumNonlinArray),
+    BSpline(Cubic(Line(OnGrid()))))
+    InterpPmm = Interpolations.extrapolate(InterpPmm, Line())
+    test_power_spectrum = 10 .^InterpPmm(log10(CosmologicalGrid.KLimberArray[1, 1]),
+    CosmologicalGrid.ZArray[1])
+    CosmoCentral.InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid,
+    BackgroundQuantitiesLoaded, PowerSpectrum, CosmoCentral.BSplineCubic())
+    @test isapprox(test_power_spectrum,
+    PowerSpectrum.InterpolatedPowerSpectrum[1, 1], rtol=1e-2)
+    CosmoCentral.WritePowerSpectrumBackground(PowerSpectrum,
+    BackgroundQuantitiesLoaded, CosmologicalGrid, "new_p_mm")
+    NewPowerSpectrum, NewBackgroundQuantitiesLoaded, NewCosmologicalGrid =
+    CosmoCentral.ReadPowerSpectrumBackground(
+    "new_p_mm",
+    MultipolesArray)
+    @test isapprox(PowerSpectrum.PowerSpectrumNonlinArray,
+    NewPowerSpectrum.PowerSpectrumNonlinArray, rtol=1e-2)
+    @test isapprox(PowerSpectrum.PowerSpectrumLinArray,
+    NewPowerSpectrum.PowerSpectrumLinArray, rtol=1e-2)
+    @test isapprox(BackgroundQuantitiesLoaded.HZArray,
+    NewBackgroundQuantitiesLoaded.HZArray, rtol=1e-2)
+    @test isapprox(BackgroundQuantitiesLoaded.rZArray,
+    NewBackgroundQuantitiesLoaded.rZArray, rtol=1e-2)
+end
+
+@testset "Test Angular coefficients evaluation" begin
+    AngularCoefficientsLoaded = CosmoCentral.ReadAngularCoefficients(
+    "/home/runner/work/CosmoCentral.jl/CosmoCentral.jl/test/cl")
+    MultipolesArray = Array(LinRange(10.5, 2999.5, 2990))
+    PowerSpectrum, BackgroundQuantities, CosmologicalGrid =
+    CosmoCentral.ReadPowerSpectrumBackground(
+    "/home/runner/work/CosmoCentral.jl/CosmoCentral.jl/test/p_mm",
+    MultipolesArray)
+    CosmoCentral.ComputeLimberArray(CosmologicalGrid, BackgroundQuantities)
+    CosmoCentral.InterpolateAndEvaluatePowerSpectrum(CosmologicalGrid,
+    BackgroundQuantities, PowerSpectrum, CosmoCentral.BSplineCubic())
+    AngularCoefficients = CosmoCentral.AngularCoefficientsStruct(
+    AngularCoefficientsArray = zeros(length(CosmologicalGrid.MultipolesArray),
+    length(GCWeightFunction.WeightFunctionArray[:, 1]),
+    length(GCWeightFunction.WeightFunctionArray[:, 1])))
+    CosmoCentral.ComputeAngularCoefficients(AngularCoefficients,
+    GCWeightFunction, GCWeightFunction, BackgroundQuantities, w0waCDMCosmology,
+    CosmologicalGrid, PowerSpectrum, CosmoCentral.CustomTrapz())
+    @test isapprox(AngularCoefficientsLoaded.AngularCoefficientsArray,
+    AngularCoefficients.AngularCoefficientsArray, rtol=1e-9)
+    CosmoCentral.WriteAngularCoefficients(AngularCoefficients, CosmologicalGrid,
+    GCWeightFunction, PiecewiseBias,
+    ConvolvedDensity, "new_cl")
+    AngularCoefficientsLoaded = CosmoCentral.ReadAngularCoefficients(
+    "new_cl")
+    @test isapprox(AngularCoefficientsLoaded.AngularCoefficientsArray,
+    AngularCoefficients.AngularCoefficientsArray, rtol=1e-9)
 end