Merge pull request #14 from JaimeRZP/cubic

Cubic Splines

src/theory.jl

@@ -1,7 +1,9 @@
 function Theory(cosmology::Cosmology,
                 names, types, pairs,
                 idx, files;
-                Nuisances=Dict())
+                Nuisances=Dict(),
+                res_wl=350, res_gc=1000,
+                int_wl="linear", int_gc="linear")
 
     nui_type =  eltype(valtype(Nuisances))
     if !(nui_type <: Float64) & (nui_type != Any)
@@ -22,7 +24,8 @@ function Theory(cosmology::Cosmology,
             zs = get(Nuisances, string(name, "_", "zs"), zs_mean)
             dz = get(Nuisances, string(name, "_", "dz"), 0.0)
             tracer = NumberCountsTracer(cosmology, zs .+ dz, nz;
-                                        b=b)
+                                        b=b, res=res_gc, 
+                                        nz_interpolation=int_gc)
         elseif t_type == "galaxy_shear"
             zs_mean, nz_mean = files[string("nz_", name)]
             m = get(Nuisances, string(name, "_", "m"), 0.0)
@@ -32,7 +35,9 @@ function Theory(cosmology::Cosmology,
             zs = get(Nuisances, string(name, "_", "zs"), zs_mean)
             dz = get(Nuisances, string(name, "_", "dz"), 0.0)
             tracer = WeakLensingTracer(cosmology, zs .+ dz, nz;
-                                       m=m, IA_params=IA_params)
+                                       m=m, IA_params=IA_params,
+                                       res=res_wl,
+                                       nz_interpolation=int_wl)
 
         elseif t_type == "cmb_convergence"
             tracer = CMBLensingTracer(cosmology)
@@ -90,7 +95,8 @@ end
 """
 function Theory(cosmology::Cosmology,
                 instructions::Instructions, files;
-                Nuisances=Dict())
+                Nuisances=Dict(),
+                kwargs...)
 
     names = instructions.names
     types = instructions.types
@@ -100,5 +106,6 @@ function Theory(cosmology::Cosmology,
     return Theory(cosmology::Cosmology,
                   names, types, pairs,
                   idx, files;
-                  Nuisances=Nuisances)
+                  Nuisances=Nuisances,
+                  kwargs...)
  end

src/tracers.jl

@@ -27,10 +27,10 @@ Kwargs:
 Returns:
 - `NumberCountsTracer::NumberCountsTracer` : Number counts tracer structure.
 """
-NumberCountsTracer(cosmo::Cosmology, z_n, nz; b=1.0, res=1000) = begin
-    nz_int = linear_interpolation(z_n, nz, extrapolation_bc=0)
-    z_w = range(0.00001, stop=z_n[end], length=res)
-    nz_w = nz_int(z_w)
+NumberCountsTracer(cosmo::Cosmology, z_n, nz; 
+    b=1.0, res=1000, nz_interpolation="linear") = begin
+    z_w, nz_w = nz_interpolate(z_n, nz, res;
+        mode=nz_interpolation)
     nz_norm = integrate(z_w, nz_w, SimpsonEven())
 
     chi = cosmo.chi(z_w)
@@ -62,16 +62,14 @@ struct WeakLensingTracer <: Tracer
 end
 
 WeakLensingTracer(cosmo::Cosmology, z_n, nz;
-    IA_params = [0.0, 0.0], m=0.0, res=350, kwargs...) = begin
-    nz_int = linear_interpolation(z_n, nz, extrapolation_bc=0) 
+    IA_params = [0.0, 0.0], m=0.0,
+    res=350, nz_interpolation="linear") = begin
     cosmo_type = cosmo.settings.cosmo_type
-    z_w = range(0.00001, stop=z_n[end], length=res)
-    dz_w = (z_w[end]-z_w[1])/res
-    nz_w = nz_int(z_w)
-    chi = cosmo.chi(z_w)
-
-    #nz_norm = sum(0.5 .* (nz_w[1:res-1] .+ nz_w[2:res]) .* dz_w)
+    z_w, nz_w = nz_interpolate(z_n, nz, res;
+        mode=nz_interpolation)
+    res = length(z_w)
     nz_norm = integrate(z_w, nz_w, SimpsonEven())
+    chi = cosmo.chi(z_w)
 
     # Calculate chis at which to precalculate the lensing kernel
     # OPT: perhaps we don't need to sample the lensing kernel
@@ -156,4 +154,24 @@ function get_IA(cosmo::Cosmology, zs, IA_params)
     A_IA = IA_params[1]
     alpha_IA = IA_params[2]
     return @. A_IA*((1 + zs)/1.62)^alpha_IA * (0.0134 * cosmo.cpar.Ωm / cosmo.Dz(zs))
-end
+end
+
+function nz_interpolate(z, nz, res; mode="linear")
+    if mode!="none"
+        if mode=="linear"
+            nz_int = linear_interpolation(z, nz;
+                extrapolation_bc=Line())
+        end
+        if mode=="cubic"
+            dz = mean(z[2:end] - z[1:end-1])
+            z_range = z[1]:dz:z[end]
+            nz_int = cubic_spline_interpolation(z_range, nz;
+                extrapolation_bc=Line())
+        end
+        zz_range = range(0.00001, stop=z[end], length=res)
+        nzz = nz_int(zz_range)
+        return zz_range, nzz
+    else
+        return z, nz
+    end
+end