diff --git a/Project.toml b/Project.toml
index 4a6b99f..98d68a1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,7 +18,7 @@ Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Lux = "b210s8857-7c20-44ae-9111-449ecde12c47"
+Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
diff --git a/data.dat b/data.dat
new file mode 100644
index 0000000..fb0f5df
Binary files /dev/null and b/data.dat differ
diff --git a/data_test.jld2 b/data_test.jld2
new file mode 100644
index 0000000..eeccae0
Binary files /dev/null and b/data_test.jld2 differ
diff --git a/examples/SIA2D.jl b/examples/SIA2D.jl
new file mode 100644
index 0000000..c8d881b
--- /dev/null
+++ b/examples/SIA2D.jl
@@ -0,0 +1,91 @@
+using OrdinaryDiffEq
+# using DifferentialEquations
+
+####  Non-allocating functions  ####
+
+nx, ny = 200, 200   # elements in the spatial grid
+Δx, Δy = 100, 100     # size in meters of each box in the grid
+n = 3.0
+ρ = 900.0
+g = 9.81
+A = 8e-18
+Γ = 2.0 * A * (ρ * g)^n / (n+2) # 1 / m^3 s 
+
+diff_x(I, Δx) = (I[begin+1:end,:] - I[1:end-1,:]) / Δx
+diff_y(I, Δy) = (I[:,begin+1:end] - I[:,1:end - 1]) / Δy
+avg(I) = (I[1:end-1,1:end-1] + I[2:end,1:end-1] + I[1:end-1,2:end] + I[2:end,2:end]) * 0.25
+avg_x(I) = (I[1:end-1,:] + I[2:end,:]) * 0.5
+avg_y(I) = (I[:,1:end-1] + I[:,2:end]) * 0.5
+inn(I) = I[2:end-1,2:end-1]
+
+year2seg = 365.25 * 24 * 60 * 60
+# tspan = [0.0, 10 * year2seg]
+tspan = [0.0, 30.0]
+
+B = zeros((nx, ny))  # elevation of the bed of the glacier
+H₀ = [ .15 * ( 70^2 - ( (i - nx/2)^2 + (j - ny/2)^2 )  ) for i in 1:nx, j in 1:ny ] # initial ice thickness
+H₀ .= max.(0.0, H₀)
+dH = zeros(size(H₀)) 
+p = nothing
+
+function SIA!(dH, H, p, t)
+    S = B .+ H
+
+    # All grid variables computed in a staggered grid
+    # Compute surface gradients on edges
+    dSdx = diff_x(S, Δx)  
+    dSdy = diff_y(S, Δy) 
+    ∇Sx = avg_y(dSdx)
+    ∇Sy = avg_x(dSdy)
+    ∇S  = (∇Sx.^2 .+ ∇Sy.^2).^((n - 1)/2) 
+
+    H̄ = avg(H)
+    
+    D = Γ .* H̄.^(n + 2) .* ∇S
+
+    # Compute flux components
+    dSdx_edges = diff_x(S[:,2:end - 1], Δx)
+    dSdy_edges = diff_y(S[2:end - 1,:], Δy)
+
+    # Cap surface elevaton differences with the upstream ice thickness to 
+    # imporse boundary condition of the SIA equation
+    # η₀ = params.physical.η₀
+    # dSdx_edges .= @views @. min(dSdx_edges,  η₀ * H[2:end, 2:end-1]/Δx)
+    # dSdx_edges .= @views @. max(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]/Δx)
+    # dSdy_edges .= @views @. min(dSdy_edges,  η₀ * H[2:end-1, 2:end]/Δy)
+    # dSdy_edges .= @views @. max(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]/Δy)
+
+    Dx = avg_y(D)
+    Dy = avg_x(D)
+    Fx = -Dx .* dSdx_edges
+    Fy = -Dy .* dSdy_edges 
+
+    #  Flux divergence
+    Fxx = diff_x(Fx, Δx)
+    Fyy = diff_y(Fy, Δy)
+    dH .= zeros(size(H))
+    dH[2:end-1, 2:end-1] .= -(Fxx .+ Fyy) 
+    @show t
+end
+
+prob = ODEProblem(SIA!, H₀, tspan, p)
+sol = solve(prob, RDPK3Sp35(), reltol = 1e-8, abstol = 1e-8, saveat=[0.0, 15.0, 30.0])
+
+using Plots
+
+H₁ = sol.u[end]
+
+
+heatmap_diff = Plots.heatmap(H₀ .- H₁, title="Variation in Ice Thickness", aspect_ratio=1, colorrange=(-200,200), colormap=:balance, clim=(-200, 200), axis=([], false))
+Plots.savefig(heatmap_diff, "variation.pdf")
+
+surface = Plots.surface(Δx * collect(1:nx), Δy * collect(1:ny), H₀, alpha=0.5, color=:green)
+
+x = Δx * collect(1:nx) / 1000 # [km]
+y1 = H₀[:, 100]
+y2 = sol.u[2][:,100]
+y3 = H₁[:, 100]
+line_plot = plot(x, [y1 y2 y3], label=["Initial condition (t=0)" "Intermediate solution (t=15)" "Final solution (t=30)"], title="Example of solutions of SIA equation", linewidth=3)
+xlabel!("Distance (km)")
+ylabel!("Heigth (m)")
+Plots.pdf(line_plot, "line_plot.pdf")
diff --git a/examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv b/examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv
new file mode 100644
index 0000000..18ff110
--- /dev/null
+++ b/examples/Torsvik_2012/Torsvik-etal-2012_dataset.csv
@@ -0,0 +1,595 @@
+Q,a95,Plate,Plate_code,Lat,Lon,CLat,CLon,RLat,RLon,Eplat,Eplong,Epang,Age
+5,4.8,north_america,101,-86.3,5.7,,,,,,,,0.5
+5,9.1,north_america,101,-86.4,8.4,,,-86.4,9.2,79.2,23,0.2,0.8
+5,7.1,north_america,101,-85.3,265.9,,,-85.3,-94.1,79.9,22.7,0.3,1
+5,8.7,north_america,101,-88.4,225.5,,,-88.4,-135.7,79.9,22.7,0.3,1
+5,4.3,north_america,101,-88.9,285,,,-88.9,-75.1,79.9,22.7,0.3,1
+3,4.9,north_america,101,-77.9,303.7,-80,3.5,-80,4.2,80.8,22.8,0.4,1.5
+3,7,north_america,101,-85.7,122.7,-81.2,73.9,-81.1,74.5,80.8,22.8,0.4,1.5
+5,5.7,north_america,101,-84.5,241,,,-84.5,-119.1,80.8,22.8,0.4,1.5
+5,6.2,north_america,101,-89.5,214.8,,,-89.5,-157.9,80.8,22.9,0.8,3
+5,4.1,north_america,101,-84.6,332.3,,,-84.8,-25.6,80.9,22.8,1,4
+5,6.7,north_america,101,-87.9,275.9,,,-88,-84.4,80.9,22.8,1,4
+5,5,north_america,101,-85.5,197.4,,,-85.5,-163.9,80.9,22.8,1.3,5
+5,5,north_america,101,-85.3,262.1,,,-85.5,-97.9,80.9,22.8,1.8,7
+5,6.9,north_america,101,-85.9,357.7,,,-86,3.9,80.9,22.8,2,8
+5,7.8,north_america,101,-82.9,309.1,,,-83.3,-47.3,80.9,22.9,2.6,10
+3,12.9,north_america,101,-81.1,225.3,-83.1,93.2,-82.4,98.8,80.9,23.2,4.2,15.5
+5,3.7,north_america,101,-82.9,346.2,,,-83.5,-1.3,80.3,25.4,5.8,21
+5,9.8,north_america,101,-82.7,301.7,,,-83.7,-51,80.2,26,6,21.7
+4,6.7,north_america,101,-87.1,189.5,,,-86.7,177.4,80,26.6,6.2,22.5
+3,8.4,north_america,101,-76.4,30.3,,,-76.2,41,79.9,26.9,6.3,23
+4,5.2,north_america,101,-80.9,331.2,,,-81.8,-17.8,79.8,27.3,6.4,23.5
+4,6.8,north_america,101,-79.7,342.8,,,-80.5,-4.3,79.3,27.5,7.1,26
+4,5.4,north_america,101,-80.2,315.4,,,-81.4,-32.9,78.8,22.7,7.5,27
+3,5,north_america,101,-82.8,316.2,,,-84.2,-24,77.4,12.5,8.6,30
+4,4.4,north_america,101,-81.9,323.6,,,-83.1,-16.9,77.4,12.5,8.6,30
+6,6.2,north_america,101,-80.6,22.1,,,-79.5,43.5,76,6.1,9.7,33
+4,4.3,north_america,101,-79.8,5.9,,,-78.9,31.9,74.9,1.4,11.4,37
+5,7,north_america,101,-78.2,297.5,,,-80.7,-40.3,74.5,0.1,12,38.5
+3,2.4,north_america,101,-85.2,296.6,,,-87.2,-9.1,74.5,-1.1,12.6,40
+6,3,north_america,101,-79.2,326.4,,,-80.1,-1.7,74.4,-2.5,13.4,42
+5,7.7,north_america,101,-83.7,323.7,-76.7,50.4,-73.2,70.3,74.3,-3.7,14.2,44
+5,8,north_america,101,-82.7,333.6,,,-82.6,19.8,74.2,-4.9,15,46
+4,9.6,north_america,101,-79.3,324.3,,,-80.2,0.8,74.2,-4.9,15,46
+4,10.1,north_america,101,-85.2,62,,,-81.2,84,74.5,-4.8,15.3,47
+6,5.6,north_america,101,-72,343.4,,,-71.9,10.7,74.5,-4.8,15.3,47
+5,4,north_america,101,-77.1,325.8,,,-78.1,-0.9,75.9,-3.5,16.2,50
+6,6.2,north_america,101,-83.1,326.3,,,-83.5,15.9,75.9,-3.5,16.2,50
+6,4.9,north_america,101,-77.6,309.1,-78.3,22.7,-75.8,51.7,76.5,-2.8,16.5,51
+6,2.6,north_america,101,-82.7,347.2,,,-81.8,31.2,76.5,-2.8,16.5,51
+6,4.4,north_america,101,-81.4,347.7,-73.9,48.8,-71.3,71.7,79.8,4.1,17.6,55
+5,14,north_america,101,-68.1,9.4,,,-67.3,34.2,80.9,6.2,18.2,57
+6,3,north_america,101,-75.9,326.4,-74.9,24.7,-73.4,52,81.8,4.8,19.4,61
+4,1.1,north_america,101,-77.1,21.3,,,-75.6,50.6,82.2,4,20,63
+4,6.6,north_america,101,-72,5.3,-68,38.1,-66.1,62.3,82.2,4,20,63
+6,3.7,north_america,101,-81.9,350.6,,,-81.7,29.2,82.2,4,20,63
+5,3.9,north_america,101,-80.5,359.4,,,-79.9,34.6,82.4,3.6,20.4,64
+5,5.8,north_america,101,-71.8,28.3,,,-69.2,56.8,81.4,-8.2,22.7,71
+5,6.2,north_america,101,-69.8,354.9,,,-68.6,26.4,81.3,-9.2,23.1,72
+6,4.6,north_america,101,-83.7,15.4,,,-80.7,58.1,80.7,-12.3,24.3,74.5
+6,9,north_america,101,-80.8,358.1,,,-78,43.4,79.5,-15.9,25.7,77
+5,9.6,north_america,101,-59.8,17.7,-57.8,41.8,-52.8,70.3,78.7,-18.1,26.9,79
+6,6.2,north_america,101,-70.5,27.6,,,-65.5,60.8,78.2,-18.8,27.5,80
+7,6.6,north_america,101,-81.8,7.8,,,-77.2,55.4,77.8,-19.4,28.1,81
+5,4.4,north_america,101,-73.6,11.1,,,-61.2,65.2,69.4,-23.5,40.5,100
+5,8.3,north_america,101,-72.4,20.7,,,-59.3,70,69.4,-23.5,40.5,100
+5,4,north_america,101,-76.4,319.4,,,-69.4,45.2,68.8,-23.1,42.6,103
+5,4.4,north_america,101,-71.5,2.4,,,-56,67.6,67.3,-22,48.2,111
+5,6.5,north_america,101,-77.4,5.9,,,-60.7,74.5,67.2,-21.8,48.9,112
+5,5.6,north_america,101,-75.7,29.1,,,-57.1,84.1,67,-21.7,49.6,113
+5,3.6,north_america,101,-74.2,30.3,,,-53.1,88.2,66,-20.6,54.2,120
+5,5.3,north_america,101,-74.8,354.8,,,-56.5,73.9,65.8,-20.2,54.8,121
+5,4.6,north_america,101,-65.9,27.8,,,-44.4,86.2,65.7,-19.8,55.3,122
+6,2.4,north_america,101,-72,11.2,,,-51.7,79.7,65.7,-19.8,55.3,122
+5,3.3,north_america,101,-71.3,7.5,,,-51.5,77.6,65.7,-19.8,55.3,122
+5,7.5,north_america,101,-71,16.9,,,-49.4,83.7,65.4,-18.9,56.8,125
+5,3.6,north_america,101,-67.2,30.8,,,-43.6,91.6,64.9,-17.7,58.7,129
+7,2.6,north_america,101,-58,23.1,,,-33.4,89.2,64.9,-16.7,63.2,144
+6,4.1,north_america,101,-64.2,338.8,-65.5,12.5,-41.3,86.4,64.8,-16.8,64.4,148
+5,3.6,north_america,101,-63.4,344.8,-64,17.5,-39.3,89.3,65.1,-16.1,64.9,150
+6,7,north_america,101,-58.7,315.1,,,-49.8,50.8,65.3,-15.7,65.2,151
+5,5.3,north_america,101,-57.4,327.5,-62.5,351,-42.6,76.4,65.8,-14.2,66.8,156
+6,6,north_america,101,-56.8,328.1,-58.1,350.6,-37.8,76.3,65.8,-13.5,69.6,163
+6,7.4,north_america,101,-52.6,318.2,-58.6,334.6,-42,69.4,65.9,-13,71.5,168
+5,4.3,north_america,101,-64,299,-71.4,321.9,-54.9,77,65.9,-13,71.5,168
+4,7.4,north_america,101,-78.7,270.3,,,-68.5,87.8,65.9,-12.9,71.9,169
+5,7.8,north_america,101,-58,303.8,,,-51.9,53,65.6,-12.9,72.9,172
+4,1.4,north_america,101,-63.2,283.2,,,-61.7,54.7,64.8,-13.4,75,180
+3,1,north_america,101,-75.7,264.7,,,-67.6,85,64.4,-13.6,75.7,183
+6,3.3,north_america,101,-58.8,256.4,-62,257.2,-72,47.8,64.1,-13.8,76.5,186
+5,7,north_america,101,-61.3,264.8,-67.8,268.3,-67.4,64.9,64.1,-13.8,76.5,186
+4,9.8,north_america,101,-58.7,277.8,-64.4,283,-61,59.4,64.1,-13.8,76.5,186
+3,7.2,north_america,101,-73.1,275.4,,,-64.1,78.6,63.8,-14,77.3,189
+5,3.1,north_america,101,-72.8,268.1,,,-66.1,79.2,63.7,-14,77.6,190
+5,8.9,north_america,101,-66,266,,,-66.7,63.8,63.3,-14.2,78.6,194
+4,2.3,north_america,101,-63,263.1,,,-67.6,56.7,63.2,-14.2,79,197
+4,11.1,north_america,101,-65.5,267.5,,,-65.9,63.2,63.2,-14.2,79,197
+6,2.8,north_america,101,-62.5,251,,,-73.2,55.4,63.2,-14.2,79,197
+5,4,north_america,101,-68,268.5,,,-65.4,69.2,63.2,-14.2,79,197
+5,6.2,north_america,101,-63.6,268.7,,,-65.2,58.6,63.2,-14.2,79,197
+5,6,north_america,101,-55.6,274.6,-59.8,273.3,-62.2,51.5,63.2,-14.2,79.1,198
+5,7.9,north_america,101,-61.5,234,,,-81,61.6,63.2,-14.2,79.2,199
+6,10.7,north_america,101,-66.4,252,,,-71.9,68.2,63.2,-14.1,79.2,200
+7,3.2,north_america,101,-66.6,268.2,,,-65.5,66.3,63.2,-14.1,79.3,201
+6,5,north_america,101,-67.8,275.8,,,-62.5,69.5,63.2,-14,79.5,204
+5,8,north_america,101,-58.5,256.9,-61.2,257.1,-70.3,51.8,63.2,-14,79.5,204
+6,10.7,north_america,101,-58.7,250.9,-62.8,251.2,-73,57.3,63.2,-14,79.5,204
+6,4.2,north_america,101,-57.8,259.3,-59.6,259.5,-68.9,47.8,63.2,-14,79.5,204
+6,6.5,north_america,101,-64.9,276.6,,,-61.9,63.3,63.2,-14,79.5,204
+6,2.5,north_america,101,-58.1,271.8,,,-62.3,48.4,63.2,-13.9,79.7,207
+6,5,north_america,101,-66.9,267.2,,,-65.8,67.6,63.2,-13.9,79.7,207
+6,4.7,north_america,101,-65.5,255.1,,,-70.7,65.7,63.2,-13.9,79.9,211
+7,5.6,north_america,101,-55.6,274.6,-59.9,273.3,-62.1,52.8,63.2,-13.9,79.9,211
+6,3.4,north_america,101,-56.6,255.9,-58.6,256.2,-70.3,45,63.2,-13.9,79.9,211
+6,4,north_america,101,-61.7,274.7,,,-61.9,57,63.2,-13.9,79.9,211
+5,3,north_america,101,-57.6,269.6,-63.3,266.9,-65.7,59,63.2,-13.9,79.9,211
+4,5.2,north_america,101,-61.4,282.2,,,-58.4,58.5,63.2,-13.9,79.9,212
+6,3.1,north_america,101,-60.1,277.1,,,-60.3,54.5,63.2,-13.9,79.9,214
+4,7,north_america,101,-60.1,271.8,,,-62.9,52.7,63.2,-13.9,79.9,215
+3,10,north_america,101,-59,267.6,,,-64.6,49,63.2,-13.9,79.9,215
+3,14,north_america,101,-56.1,276,-62.4,280.5,-59.5,59.9,63.2,-13.9,79.9,215
+7,7.8,north_america,101,-50.5,267.6,-53.4,268.7,-62.1,37.9,63.2,-13.9,79.9,215
+6,5.6,north_america,101,-57.4,267.7,-59.3,268.4,-64.3,50,63.2,-13.9,79.9,216
+6,3.2,north_america,101,-59.9,279.4,,,-59.2,54.9,63.2,-13.9,79.9,217
+6,5.1,north_america,101,-52.9,282,-53.5,282.5,-54.9,45.6,63.2,-13.9,79.9,218
+6,7.7,north_america,101,-56.4,276.8,-58,277.6,-59.4,50.7,63.2,-13.9,79.9,218
+7,5,north_america,101,-59.6,277.5,-63.5,280.6,-59.7,61.9,63.2,-13.9,79.9,220
+6,5,north_america,101,-53.4,281.7,-55.9,281.5,-56.6,48.9,63.2,-13.9,79.9,220
+6,2,north_america,101,-58.5,279.5,,,-58.6,52.4,63.2,-13.9,79.9,221
+5,5,north_america,101,-60.5,281.6,,,-58.4,56.7,63.2,-13.9,79.9,221
+4,3.9,north_america,101,-48.3,272.3,,,-57.6,31.5,63.2,-13.9,79.9,221
+6,5,north_america,101,-49.1,285.1,-52.7,287.8,-51.7,47.3,63.2,-13.9,79.9,225
+6,4,north_america,101,-54.1,285.2,,,-53.8,47.9,63.2,-13.9,79.9,227
+6,3.2,north_america,101,-48.4,278.5,,,-54.2,35.8,63.2,-13.9,79.9,228
+5,2.5,north_america,101,-54.1,288.3,-58.6,293.1,-52.3,58,63.2,-13.9,79.9,230
+6,4,north_america,101,-45.2,295.4,-49.1,299,-44,49.4,63.2,-13.9,79.9,230
+5,3.1,north_america,101,-52.5,290.7,-56.6,295,-50.4,56.1,63.2,-13.9,79.9,230
+6,4.5,north_america,101,-56.5,283.2,-59.8,286.6,-55.8,57.3,63.2,-13.9,79.9,230
+6,3.3,north_america,101,-46.1,293.6,-49.8,296.9,-45.5,48.9,63.2,-13.9,79.9,230
+7,3.4,north_america,101,-55.5,287.9,-58.2,290.7,-53.2,56.5,63.2,-13.9,79.9,234
+7,5,north_america,101,-54.6,284.5,-57.8,287.8,-54.4,54.5,63.2,-13.9,79.9,234
+6,4.9,north_america,101,-44.7,301.4,-44.9,301.7,-39.8,46.9,63.2,-13.9,79.9,234
+5,4.9,north_america,101,-55.6,285.8,-60.5,290.9,-54.2,59.9,63.2,-13.9,79.9,234
+4,5.3,north_america,101,-41.1,305.6,-40.8,305.2,-35,45.7,63.2,-13.9,79.9,234
+5,7,north_america,101,-46.1,301,-50.6,306.1,-41.5,55.1,63.2,-13.9,79.9,235
+6,7.2,north_america,101,-44.3,271.6,-45.5,271.1,-56.6,26.5,63.2,-13.9,79.9,246
+6,5,north_america,101,-51,306.5,-53.8,310.3,-41.8,60.5,63.2,-13.9,79.9,250
+4,8,north_america,101,-49.9,298.1,-51.6,300,-45.1,52.6,63.2,-13.9,79.9,255
+4,15,north_america,101,-54.8,299.3,-57.3,302.2,-47.6,60.4,63.2,-13.9,79.9,258
+5,5,north_america,101,-51.5,306.7,-54.8,311.6,-41.9,62.2,63.2,-13.9,79.9,263
+5,3.8,north_america,101,-56.3,302.9,,,-46.6,59.5,63.2,-13.9,79.9,270
+4,3.8,north_america,101,-53,308.7,,,-42,58.9,63.2,-13.9,79.9,272
+4,8.6,north_america,101,-54.8,292.1,,,-50.8,52.3,63.2,-13.9,79.9,272
+3,10,north_america,101,-51.9,303,-56.5,311,-43.2,63.6,63.2,-13.9,79.9,277
+3,5,north_america,101,-51.7,302.1,-53.7,304,-44.5,57.1,63.2,-13.9,79.9,277
+7,3.6,north_america,101,-42.1,306.5,-41.4,306.4,-34.9,46.9,63.2,-13.9,79.9,280
+5,16.3,north_america,101,-33.5,306.3,,,-29,40.8,63.2,-13.9,79.9,282
+4,13.1,north_america,101,-44.6,305.3,-47.8,308.8,-38.3,54,63.2,-13.9,79.9,283
+5,2.1,north_america,101,-46.8,304,-48.3,305.9,-40.1,52.7,63.2,-13.9,79.9,285
+3,10.2,north_america,101,-38.9,300.8,,,-35.9,41,63.2,-13.9,79.9,289
+5,1.5,north_america,101,-50.5,303,-56.2,310.3,-43.4,63,63.2,-13.9,79.9,291
+4,5,north_america,101,-37.5,296.6,-35.8,295.2,-36.7,34.4,63.2,-13.9,79.9,292
+5,7.1,north_america,101,-40.1,307.7,-42.1,310.1,-33.5,50,63.2,-13.9,79.9,292
+5,2,north_america,101,-43.1,307.9,-46.5,311.7,-36,54.6,63.2,-13.9,79.9,292
+4,2,north_america,101,-41.6,300.4,-42.6,301.5,-38.3,44.7,63.2,-13.9,79.9,292
+4,2.8,north_america,101,-40.1,300.5,-40.3,300.8,-37,42.2,63.2,-13.9,79.9,298
+4,12.8,north_america,101,-55.3,279.8,-60.6,285.1,-56.8,58,63.2,-13.9,79.9,299
+5,3.9,north_america,101,-44.1,301.5,-41.5,300.4,-38,43,63.2,-13.9,79.9,300
+5,2.1,north_america,101,-42.1,312.1,-43,313.4,-32.7,52.9,63.2,-13.9,79.9,301
+5,3.4,north_america,101,-44.1,303.9,-46.3,306.8,-38.2,51.4,63.2,-13.9,79.9,301
+6,3.1,north_america,101,-28.6,299.9,,,-28.6,32.4,63.2,-13.9,79.9,303
+5,1.8,north_america,101,-45.7,308.6,-50.5,314.6,-37.6,59.8,63.2,-13.9,79.9,303
+5,6,north_america,101,-36,302,-30.2,301.5,-29,34.8,63.2,-13.9,79.9,310
+7,7.7,north_america,101,-27.2,298.3,,,-28.4,30.2,63.2,-13.9,79.9,317
+6,8.3,north_america,101,-22.6,294.4,,,-26.9,23.8,63.2,-13.9,79.9,320
+7,15.3,north_america,101,-27.9,297.2,,,-29.5,29.8,63.2,-13.9,79.9,322
+4,6.5,north_america,101,-19.5,315.8,,,-12.6,39.2,63.2,-13.9,79.9,330
+6,8,north_america,101,-27,311,-17.8,309.8,,,,,,333
+7,9,north_america,101,-18.6,304.2,,,,,,,,335
+3,16,north_america,101,-27.4,303,-16.6,299.6,,,,,,370
+4,9,north_america,101,-13,285,1.5,284.8,,,,,,415
+7,5.3,north_america,101,-17,305,,,,,,,,420
+6,5.8,north_america,101,-19.1,308.3,,,,,,,,425
+4,7.3,north_america,101,-24,326.6,-16.9,321.7,,,,,,438
+4,3.9,north_america,101,-13.4,329.3,,,,,,,,470
+4,4.3,north_america,101,-17.5,332.4,,,,,,,,480
+6,11.9,north_america,101,-10.3,346.5,,,,,,,,490
+5,8.5,north_america,101,0.6,343,3.1,338.9,,,,,,495
+5,9.7,north_america,101,-10.6,338,-8.4,334.6,,,,,,495
+5,9,north_america,101,-5.2,345.8,-4.7,345,,,,,,495
+5,7.1,north_america,101,3.4,355.1,,,,,,,,500
+6,4.3,north_america,101,-12.6,337.3,,,,,,,,500
+6,10,north_america,101,5.4,348.7,,,,,,,,503
+5,3.3,north_america,101,-0.6,341.1,-1.7,342.6,,,,,,508
+5,6.2,north_america,101,11.9,4.5,,,,,,,,532
+6,7.4,greenland,102,-76.3,21.5,,,-75,44.1,72.8,9.1,11.2,39
+5,6.3,greenland,102,-74.6,339.4,,,-77.6,11.6,71.1,30.9,16,54
+4,6,greenland,102,-63,0,,,-64.6,25.6,71.4,30.2,16.4,54.5
+4,15,greenland,102,-63,354,,,-65,19.8,71.7,29.6,16.8,55
+3,4.2,greenland,102,-61,345,,,-63.8,9.9,72,29,17.2,55.5
+4,8.9,greenland,102,-63.4,5.1,,,-64.4,32,72,29,17.2,55.5
+6,3,greenland,102,-68,358,,,-69.6,28,72,30.1,17.5,59
+5,9,greenland,102,-56,3,,,-57.4,27.7,71.9,30.3,17.5,59.5
+6,3.2,greenland,102,-67.5,15,,,-67.6,44.9,71.8,30.8,17.5,60.5
+5,6,greenland,102,-76.2,37.9,,,-73.9,73.2,71.8,30.8,17.5,60.5
+5,6.9,greenland,102,-64.8,321.5,,,-69.4,-14.9,71.8,30.8,17.5,60.5
+6,5,greenland,102,-52.7,278.7,,,-67.9,36.6,60.5,1.4,69.3,209
+3,11,greenland,102,1,302,,,,,,,,400
+4,3.6,europe,301,-80.6,267.5,,,-80.6,-92.3,18.3,-47,0.1,0.5
+4,4.4,europe,301,-86.4,296.1,,,-86.4,-63.1,18.3,-47,0.1,0.5
+3,12.9,europe,301,-84.3,357.7,,,-83.9,5.8,18,-26.7,1,8
+5,1.8,europe,301,-78.9,328.3,,,-78.9,-25.8,17.9,-26.4,1.1,9.5
+5,3.5,europe,301,-77.4,314.2,,,-77.7,-40.4,18.5,-26.3,1.2,10
+5,6.9,europe,301,-84.1,351.2,,,-83.6,3,19.7,-25.7,1.4,11.5
+5,1.5,europe,301,-72.4,352,,,-71.9,-3.5,19.6,-25.5,1.5,12
+4,4.4,europe,301,-77.8,310.8,,,-78.9,-36.2,20,-19.7,2.9,24
+3,3.4,europe,301,-80.8,2,,,-78,26.5,26.8,-19.5,5.5,34
+5,6.8,europe,301,-63.7,358.6,-76.7,357.2,-72.2,30.6,30.5,-15.7,10.3,52
+5,1.5,europe,301,-83,335,,,-78.2,41.3,34.5,-15.7,12.4,59
+5,3.5,europe,301,-76,340,,,-72.9,22.7,34.5,-15.7,12.4,59
+5,1.2,europe,301,-81.7,359.8,,,-74.9,46.3,34.6,-15.7,12.5,59.4
+5,2.8,europe,301,-80.2,339.6,,,-76,33.2,34.6,-15.7,12.5,59.5
+4,1.5,europe,301,-82.5,338,,,-77.6,40.6,34.7,-15.7,12.5,59.8
+5,5.7,europe,301,-77.7,325.4,,,-76.2,18,34.7,-15.7,12.5,59.8
+4,2.7,europe,301,-77,355,,,-71.5,34.4,34.8,-15.7,12.6,60
+6,4.5,europe,301,-71.4,334.7,,,-69.8,11.3,34.8,-15.7,12.6,60.1
+6,2.4,europe,301,-81,359,,,-74.3,44.8,34.9,-15.7,12.7,60.7
+6,4.7,europe,301,-73.3,346.2,,,-69.5,23.1,35,-15.7,12.8,60.9
+5,4.7,europe,301,-78.9,347,,,-74,34.1,35,-15.7,12.8,61
+5,2.7,europe,301,-74,351,,,-69.4,27.5,35.1,-15.7,12.8,61.2
+6,3.4,europe,301,-79.1,344.5,-84.9,222.9,-82.4,95.9,35.8,-16,14.6,68
+7,8,europe,301,-73,336,,,-69.8,21.1,35.4,-15.7,15.2,74
+5,3,europe,301,-74,341,,,-67.7,35.4,37.5,-14,19.5,86
+5,3,europe,301,-74,328,,,-70.4,28.7,37.5,-14,19.5,86
+4,5,europe,301,-68,329,,,-65.1,22.3,39.2,-13.6,21.8,89.5
+5,4,europe,301,-76,1,,,-62.9,54.7,39.4,-13.7,24.2,93
+6,2.5,europe,301,-80.8,338.4,,,-63.2,70.6,41.6,-10.1,34.1,108
+6,2.9,europe,301,-74,3,,,-48,82.4,45.9,-3.4,48.4,140
+5,6,europe,301,-78,328,,,-56,87.9,48.4,-1.1,52.8,156.5
+3,3.9,europe,301,-70,327,,,-54.5,74.3,48.5,-1,53,157
+5,7,europe,301,-78,310,,,-59.5,88.7,48.6,-0.9,53.4,158
+4,4,europe,301,-72,312,,,-59.3,77.2,48.7,-0.8,53.8,159
+6,7.3,europe,301,-72,330,,,-53.9,78.9,48.7,-0.8,53.8,160
+7,6.3,europe,301,-63,300,,,-62.2,62.5,49.8,0.3,57,168
+6,6.8,europe,301,-69,283,,,-65.8,82.6,49.4,0.2,60.3,179
+6,12,europe,301,-71,276,,,-65.6,90.9,49.2,0,61.8,184
+5,5.1,europe,301,-66,295.2,,,-60.6,75.2,49.1,-0.1,62.3,186
+4,3,europe,301,-77,315,,,-52.1,94.8,48.7,0,63.6,192
+5,7.5,europe,301,-61,259,,,-76,88.7,48.1,0.7,63.1,198
+6,9,europe,301,-55,280,,,-69.1,54,47.8,1.1,62.5,201
+7,3,europe,301,-51,285,,,-65.8,44.1,47.8,1.1,62.5,201
+6,4.5,europe,301,-50,286.4,,,-65.3,41.7,47.6,1.6,61.9,204.2
+6,8,europe,301,-50,292,-58,272.9,-73.6,65.6,47.3,2.2,61.1,208
+5,5.1,europe,301,-50,308,-58.2,298.2,-61.7,60.7,46.6,3.2,59.5,215
+5,4.6,europe,301,-50,305,,,-56.9,45.8,46,3.9,58.2,221
+6,6,europe,301,-49,311,,,-52.9,46.9,46,3.9,58.2,226
+5,2.9,europe,301,-47.1,301.6,,,-57.8,39.2,46,3.9,58.2,228
+6,3,europe,301,-54,321,,,-49.4,58.3,46,3.9,58.2,234
+6,12,europe,301,-53,303,,,-59.1,50.4,46,3.9,58.2,234
+5,15,europe,301,-49,326,-56.5,318.6,-51.6,61.2,46,3.9,58.2,239
+6,5,europe,301,-43,326,-47.8,322.2,-45.8,50.6,46,3.9,58.2,243
+6,3.8,europe,301,-49,348.2,-57,344.7,-39.2,71.2,46,3.9,58.2,246
+6,7.8,europe,301,-59.3,325.8,,,-49,67.4,46,3.9,58.2,248
+7,3.3,europe,301,-50.6,345.6,-58.8,341,-41.8,71.7,46,3.9,58.2,249
+3,10,europe,301,-59,330,,,-46.9,68.2,46,3.9,58.2,250
+6,3.3,europe,301,-56.2,326,,,-47.7,63.2,46,3.9,58.2,251
+3,13.9,europe,301,-52.7,328.4,,,-44.8,59.6,46,3.9,58.2,251
+5,2.1,europe,301,-56.4,321.7,,,-49.9,62,46,3.9,58.2,251
+6,5,europe,301,-50,343,-59.3,336.5,-44.1,70.7,46,3.9,58.2,251
+3,5.3,europe,301,-53.3,330.2,,,-44.2,61.1,46,3.9,58.2,251
+3,5,europe,301,-54.3,323,,,-48.4,59.5,46,3.9,58.2,251
+3,2.7,europe,301,-58.5,314.5,,,-54.3,63.2,46,3.9,58.2,251
+6,9.7,europe,301,-52.8,334.4,,,-41.8,62.3,46,3.9,58.2,251
+5,2.7,europe,301,-46,327,-50,324.7,-45.5,54.7,46,3.9,58.2,255
+6,4,europe,301,-51,341,-55.4,338.1,-41.4,66.8,46,3.9,58.2,255
+5,3.5,europe,301,-45.6,350.2,-54.9,338.2,-41,66.2,46,3.9,58.2,260
+5,5,europe,301,-47,331,-50.6,327.7,-44.2,56.7,46,3.9,58.2,261
+6,4,europe,301,-49,343,-52.6,341.7,-38,65.3,46,3.9,58.2,261
+4,0,europe,301,-53,331,-58.6,326.4,-48.5,66.5,46,3.9,58.2,264
+5,1.5,europe,301,-49,334,-52.5,331.5,-43.1,60.7,46,3.9,58.2,264
+3,4.6,europe,301,-47,336,-49.3,334.7,-39.6,58.5,46,3.9,58.2,264
+5,4,europe,301,-51,324,-56.4,318.8,-51.4,61,46,3.9,58.2,264
+5,6.1,europe,301,-51.5,322,,,-47.7,55.3,46,3.9,58.2,264
+4,4.1,europe,301,-50,344,,,-35.3,63.8,46,3.9,58.2,269
+5,2.5,europe,301,-51,343,,,-36.4,64.3,46,3.9,58.2,271
+5,5.9,europe,301,-53,344,,,-37.1,66.7,46,3.9,58.2,271
+4,8.6,europe,301,-51,345,,,-35.4,65.3,46,3.9,58.2,275
+4,11,europe,301,-54,352,,,-34.1,71.4,46,3.9,58.2,279
+4,7,europe,301,-37,341,,,-28.3,50.6,46,3.9,58.2,280
+3,14,europe,301,-47,337,,,-37.1,57.3,46,3.9,58.2,280
+5,10,europe,301,-42,346,,,-29,57.8,46,3.9,58.2,280
+5,1,europe,301,-47,337,,,-37.1,57.3,46,3.9,58.2,281
+4,13.4,europe,301,-44.6,337.4,,,-35.4,55.2,46,3.9,58.2,281
+4,6.9,europe,301,-38,346,,,-26.2,54.6,46,3.9,58.2,281
+5,6.5,europe,301,-49.4,359.7,,,-27.4,71.4,46,3.9,58.2,282.6
+4,6.7,europe,301,-41,352,,,-25.1,60.6,46,3.9,58.2,285
+5,3.2,europe,301,-43,352,,,-26.5,62.1,46,3.9,58.2,285
+5,5.1,europe,301,-44,4,-48.6,4.1,-24.9,73.2,46,3.9,58.2,285
+5,2,europe,301,-40,346,-42.4,345.2,-29.6,57.7,46,3.9,58.2,285
+3,7.7,europe,301,-44,350,-39.8,-9.4,-25,58.8,46,3.9,58.2,285
+5,6.3,europe,301,-38,346,,,-26.2,54.6,46,3.9,58.2,285
+4,8.1,europe,301,-42,354,,,-24.8,62.5,46,3.9,58.2,285
+5,2,europe,301,-42,349,-40.5,-9.8,-25.7,59,46,3.9,58.2,285
+4,17,europe,301,-49,342,-52.5,340.1,-38.7,64.5,46,3.9,58.2,285
+4,6.8,europe,301,-37,340,-38.5,339.4,-30.2,50.9,46,3.9,58.2,285
+4,7.9,europe,301,-43,345,,,-30.2,58.1,46,3.9,58.2,285
+5,4,europe,301,-41,345,-44,343.5,-31.6,58.1,46,3.9,58.2,285
+3,13.2,europe,301,-40,352,,,-24.3,59.8,46,3.9,58.2,285
+5,4,europe,301,-50,330,,,-42.6,57.1,46,3.9,58.2,286
+4,5.9,europe,301,-49,356,,,-28.8,69.1,46,3.9,58.2,286
+4,10,europe,301,-48,343,,,-34.5,61.4,46,3.9,58.2,286
+3,1,europe,301,-42,353,,,-25.3,61.9,46,3.9,58.2,286
+4,5.8,europe,301,-41.5,340,-45.3,338.1,-35.5,56.3,46,3.9,58.2,287
+4,2.4,europe,301,-32,354,,,-17.3,55.4,46,3.9,58.2,291
+3,15.9,europe,301,-41,349,,,-26.7,58.8,46,3.9,58.2,291
+3,13,europe,301,-46,347,,,-31.2,61.8,46,3.9,58.2,291
+6,13,europe,301,-42,346,,,-29,57.8,46,3.9,58.2,293
+4,4.8,europe,301,-44,339,,,-34.1,55.5,46,3.9,58.2,294
+5,3.5,europe,301,-32.9,347.1,,,-21.9,51.4,46,3.9,58.2,294
+5,6.3,europe,301,-35.4,346.8,,,-23.9,53.1,46,3.9,58.2,294
+4,4,europe,301,-47,348,,,-31.3,63.2,46,3.9,58.2,294
+3,19,europe,301,-42,348,,,-27.9,59,46,3.9,58.2,294
+4,4.8,europe,301,-44,355,,,-25.7,64.6,46,3.9,58.2,294
+5,8.1,europe,301,-47.1,337.1,,,-37.1,57.5,46,3.9,58.2,294
+5,6.5,europe,301,-38,348,,,-25.1,55.8,46,3.9,58.2,294
+5,11,europe,301,-37,354,,,-21.1,58.9,46,3.9,58.2,294
+4,7.1,europe,301,-37.1,350,,,-23.3,56.4,46,3.9,58.2,295
+4,13.6,europe,301,-43,354,,,-25.5,63.3,46,3.9,58.2,296
+3,7.1,europe,301,-42.5,339.6,,,-32.8,54.5,46,3.9,58.2,297
+4,2.9,europe,301,-39,341,,,-29.7,52.3,46,3.9,58.2,297
+5,1.3,europe,301,-41,342,,,-30.5,54.6,46,3.9,58.2,297
+6,3,europe,301,-43,345,-49.9,337.3,-38.6,60.4,46,3.9,58.2,297
+5,2.4,europe,301,-48.4,349.8,-56.1,341.2,-40.3,68.8,46,3.9,58.2,299
+5,4,europe,301,-31,354,,,-16.6,54.7,46,3.9,58.2,299
+5,2.2,europe,301,-48.2,348.3,-55.9,339.4,-41,67.9,46,3.9,58.2,299
+6,4,europe,301,-42,359,-45.7,356.5,-26.3,66.8,46,3.9,58.2,301
+5,2,europe,301,-48.2,342.3,,,-35,61.3,46,3.9,58.2,303
+3,3,europe,301,-49,349,,,-32.2,65.4,46,3.9,58.2,303
+5,5.2,europe,301,-38.3,354,,,-22,59.8,46,3.9,58.2,305
+5,9,europe,301,-38,343,-40.3,341.5,-30.3,53.7,46,3.9,58.2,305
+5,2.9,europe,301,-38.4,339.5,,,-30.1,50.9,46,3.9,58.2,312
+6,7,europe,301,-14,332,,,,,,,,332
+6,8.2,europe,301,-14.3,335.9,,,,,,,,335
+6,11,europe,301,4,323,,,,,,,,396
+6,7,europe,301,-5,320,,,,,,,,410
+7,3,europe,301,2,318,,,,,,,,410
+6,7.3,europe,301,3.7,325.5,,,,,,,,411
+5,6,europe,301,2,321,,,,,,,,412
+4,2.5,europe,301,-8,335,,,,,,,,415
+4,4.3,europe,301,-1,325,,,,,,,,415
+4,4.9,europe,301,-17,350,,,,,,,,419
+7,9.1,europe,301,-19,344,-13.6,-16,,,,,,421
+3,5.7,europe,301,-27,346,,,,,,,,421
+3,8,europe,301,-23,351,,,,,,,,422
+4,4.6,europe,301,-15,350,,,,,,,,424
+5,7.9,europe,301,-15,347,,,,,,,,425
+3,6,europe,301,-21,344,,,,,,,,425
+4,2,europe,301,-19,349,,,,,,,,426
+5,5.1,europe,301,-19,352,,,,,,,,427
+3,5,europe,301,-21,344,,,,,,,,430
+4,5.2,europe,301,-21,358,,,,,,,,432
+5,13.4,europe,301,3,35,,,,,,,,458
+6,4.8,europe,301,5,34,,,,,,,,459
+5,4.9,europe,301,12,41.9,,,,,,,,463
+6,4.4,europe,301,14,49,,,,,,,,466
+4,6.7,europe,301,11.3,39.1,,,,,,,,467
+6,9.2,europe,301,19,51,,,,,,,,471
+6,6.8,europe,301,18.7,54,,,,,,,,472
+4,7.1,europe,301,25,50.9,,,,,,,,472
+5,9,europe,301,30,55,,,,,,,,475
+6,5.1,europe,301,18,46,,,,,,,,475
+6,4,europe,301,18,55,,,,,,,,477
+5,3.6,europe,301,34.7,59.1,,,,,,,,478
+5,5,europe,301,22,87,34.2,79.8,,,,,,500
+3,6.8,europe,301,52,111,,,,,,,,500.1
+4,6.9,europe,301,58.4,122.5,68.7,102.2,,,,,,535
+3,10,Amazonia,201,-85.4,303.8,,,-85.5,-37.4,62.1,-40.5,2.8,8.5
+3,10,Amazonia,201,-85.7,80.5,,,-84.6,75.6,62.1,-40.5,2.8,8.5
+5,11.4,Amazonia,201,-84.5,316,,,-84.4,-26.3,62,-40.6,3.1,9.5
+5,5.9,Amazonia,201,-79.5,0,,,-69.6,46.8,63.5,-33.4,26.2,70.5
+6,2.6,Amazonia,201,-83.2,320.1,,,-71.8,48.8,61.6,-34.3,33.8,84
+5,4.2,Amazonia,201,-79.4,331.9,,,-67.5,44.3,61.2,-34.3,34.4,85
+3,4.8,Amazonia,201,-87.6,315.1,,,-69.5,65.2,58.2,-34.6,38.7,92
+4,3,Amazonia,201,-87.9,335.9,,,-63.5,69.8,54.9,-34.8,44.9,102
+4,2.8,Amazonia,201,-83.6,81,,,-53.7,84.1,51.8,-35,52.4,118
+6,2.6,Amazonia,201,-89.1,3.3,,,-56.3,76.3,51,-34.3,53.6,123.5
+6,2,Amazonia,201,-82.4,30.3,,,-48.1,77.9,50.1,-32.8,54.8,130
+6,2.4,Amazonia,201,-83,71.4,,,-49.6,85.6,50,-32.5,55.1,132
+4,14.1,Amazonia,201,-80.6,275,,,-59.3,63.3,50,-32.5,55.1,146
+3,9.3,Amazonia,201,-85.3,262.5,,,-58.5,72.7,50,-32.5,55.1,175
+5,3.8,Amazonia,201,-65.5,250.3,,,-70.6,34.3,50,-32.5,55.1,196.6
+5,4,Amazonia,201,-81.2,235.1,,,-63.6,72.7,50,-32.5,55.1,198
+4,4.9,Amazonia,201,-66.9,245.6,,,-71.9,40.3,50,-32.5,55.1,202.5
+4,10,Amazonia,201,-82,320,,,-52.7,66.4,50,-32.5,55.1,232
+6,6,Amazonia,201,-71.4,303.6,-60,294.1,-50.4,29.1,50,-32.5,55.1,248.5
+6,6.6,Amazonia,201,-80.7,7,-70.7,325.5,-45.5,52.9,50,-32.5,55.1,260
+4,4.5,Amazonia,201,-68.2,321.3,-56.1,305.2,-43.4,29.2,50,-32.5,55.1,280
+4,4.1,Amazonia,201,-65.7,330.9,-53.2,324,-33.7,36,50,-32.5,55.1,300
+3,11.2,Amazonia,201,-57,357,-48.3,338.2,-24.3,41.3,50,-32.5,55.1,310
+3,6,Parana,202,-77,18,,,-66.8,52.3,63.9,-33.6,24.7,65.5
+4,3.7,Parana,202,-84.6,115.4,,,-57,86.1,51.7,-35,52.6,119
+3,10.4,Parana,202,-81,14,,,-49.2,71.8,51.5,-34.9,53,121
+5,5.9,Parana,202,-72,25,,,-39.1,73.6,50.9,-34.2,53.7,124
+3,11,Parana,202,-79,8,,,-43,71.3,47.5,-33.3,56.1,139.5
+3,18.1,Parana,202,5.9,338.1,,,24.1,12.1,47.6,-33.3,56.2,510
+4,8,Colorado,290,-85,222,,,-73,71.4,56.9,-34.7,40.9,95.5
+5,8.7,Colorado,290,-83,138,,,-53.3,90.1,47.5,-33.3,57.3,183
+6,4.5,Colorado,290,-74,67,,,-38.1,89.4,47.5,-33.3,57.3,183
+5,6.8,Colorado,290,-75.5,129.4,,,-51,102,47.5,-33.3,57.3,187
+6,4.5,Colorado,290,-51,223,,,-84.7,-21.3,47.5,-33.3,57.3,195
+4,7.6,Colorado,290,-81.8,298.3,,,-52.4,65.2,47.5,-33.3,57.3,216
+4,7,Colorado,290,-83,317,-69.1,298.5,-49.7,45.5,47.5,-33.3,57.3,240
+6,6.4,Colorado,290,-76,313.4,,,-48.1,57.7,47.5,-33.3,57.3,240
+6,4.9,Colorado,290,-89.2,346.1,-75.1,293.5,-52.6,54.2,47.5,-33.3,57.3,240
+6,3.3,Colorado,290,-80.1,348.6,,,-44.9,68.6,47.5,-33.3,57.3,263
+6,4.1,Colorado,290,-75.7,326,,,-45.3,59.5,47.5,-33.3,57.3,264
+4,5.2,Colorado,290,-80.6,308.3,,,-50.7,63.7,47.5,-33.3,57.3,267
+4,2.8,Colorado,290,-80.6,268.8,-66.7,285.5,-53.7,39.6,47.5,-33.3,57.3,283
+3,2.5,Colorado,290,-59.5,357.5,-55.2,332.3,-29.2,43.2,47.5,-33.3,57.3,283
+3,3.1,Colorado,290,-74,313,-60.9,301.3,-45.1,35.2,47.5,-33.3,57.3,283
+4,4.9,Colorado,290,-76.8,293.7,-62.4,292.2,-49.6,34.4,47.5,-33.3,57.3,290
+4,7,Colorado,290,-73.1,272.4,,,-58.4,50.8,47.5,-33.3,57.3,290
+3,5,Colorado,290,-66,348,,,-33.4,57.8,47.5,-33.3,57.3,300
+4,5.7,Colorado,290,-51,347,,,-20.8,48.8,47.5,-33.3,57.3,310
+6,5,Colorado,290,-49,343,-45.4,323.6,-25.1,31.5,47.5,-33.3,57.3,310
+6,9.6,Colorado,290,-57,350,,,-25.2,53.6,47.5,-33.3,57.3,310
+5,5.7,Patagonia,291,-81,337.4,,,-74.4,33.2,57.2,-31.2,19.4,47
+6,4.3,Patagonia,291,-78.7,358.4,,,-68.8,45.5,63.4,-33.4,26.5,71.5
+5,2,Patagonia,291,-86.8,35.2,,,-66.7,71.5,58.2,-34.6,38.7,92
+4,3.8,Patagonia,291,-87,159,-77.5,283.2,-60.9,51.4,52,-35,52.1,116
+5,5.5,Patagonia,291,-81,172,,,-57.9,90.4,47.5,-33.3,58,156.5
+4,4.9,Patagonia,291,-85,197,,,-55.8,82,47.5,-33.3,59.2,167
+5,5.2,Southern_Africa,701,-64.1,46.1,,,,,,,,90.5
+6,9.7,Southern_Africa,701,-47.6,89.9,,,,,,,,129
+5,3.1,Southern_Africa,701,-48.3,86.6,,,,,,,,132
+4,13.3,Southern_Africa,701,-64,80.6,,,,,,,,180
+3,15.8,Southern_Africa,701,-61.9,71.9,,,,,,,,183
+5,3.2,Southern_Africa,701,-71.6,93.5,,,,,,,,183
+4,11,Southern_Africa,701,-70.5,88.7,,,,,,,,183
+3,7,Southern_Africa,701,-57,84,,,,,,,,183
+5,9.5,Southern_Africa,701,-65.4,75.1,,,,,,,,183
+5,8.7,Southern_Africa,701,-70.7,106.7,,,,,,,,186
+3,4.6,Southern_Africa,701,-68,50.5,-54.7,39.5,,,,,,221.5
+3,6,Southern_Africa,701,-54,80,-49,62.6,,,,,,248.5
+5,7.6,Southern_Africa,701,-50.9,86.3,-48,63.8,,,,,,251
+6,5.6,Southern_Africa,701,-26.5,26.5,-21.6,27.7,,,,,,281.5
+6,12,Southern_Africa,701,-25,67,-25.2,53.6,,,,,,315
+4,7,Southern_Africa,701,10,15,-3.3,16,,,,,,398.5
+5,18,Southern_Africa,701,25,343,17,-11,,,,,,446.5
+6,9,Southern_Africa,701,28,14,13.9,13.8,,,,,,482.5
+4,14,Meseta,707,-46,78,,,,,,,,120
+6,9,Meseta,707,-44,71,,,-42.6,73.6,33.6,26,2.3,173.5
+4,11,Meseta,707,-45,68,,,-43.7,70.7,33.6,26,2.3,173.5
+5,6,Meseta,707,-69.2,55.5,,,-68.2,61,33.6,26,2.3,201
+3,7,Meseta,707,-71,36,,,-70.5,42.7,33.6,26,2.3,201
+5,19,Meseta,707,-73,61.3,,,-71.8,67.4,33.6,26,2.3,201
+3,4.6,Meseta,707,-38.7,56.8,,,-37.7,59.4,33.6,26,2.3,273
+5,4.7,Meseta,707,-32.2,64.1,-33.4,66,-32.1,68.3,33.6,26,2.3,273
+3,7.8,Meseta,707,-24,63.8,-23.3,62.3,-22.1,64.3,33.6,26,2.3,273
+3,20.9,Meseta,707,-36,58,,,-34.9,60.5,33.6,26,2.3,280.5
+3,4.1,Somalia,709,-87.5,359.3,,,-87.5,-2.3,52.4,6.3,-0.1,1
+3,5.7,Somalia,709,-87.2,37.1,,,-87.3,34.3,50.2,6.5,-0.2,2
+3,4.1,Somalia,709,-79.7,350.2,,,-79.6,-11,50.3,6.4,-0.3,2.5
+5,3.8,Somalia,709,-85.7,75.8,-84.4,65.5,-85.1,59.6,50.4,6.4,-1.3,11.5
+4,8.8,Somalia,709,-86.5,6.6,,,-86.4,-7.5,50.4,6.4,-1.3,12
+3,10,Somalia,709,-80.1,214.2,,,-79.7,-142.7,50.3,6.4,-1.3,13.5
+5,3.1,Somalia,709,-84.6,343.3,,,-84.2,-25.5,50.4,6.4,-1.3,17
+4,8.4,Somalia,709,-75.1,350.3,,,-74.8,-14.2,50.3,6.4,-1.5,34
+4,3,Somalia,709,-61.8,79.5,,,-62.7,77.8,50.3,6.4,-1.5,111
+3,9.3,Somalia,709,-60,82,,,-60.9,80.4,50.3,6.4,-1.5,124.5
+3,5,Somalia,709,-46,40,-33.3,36.9,-33.7,35.2,50.3,6.4,-1.5,257
+5,1.9,Somalia,709,27.5,344.8,,,27.8,-15.9,50.3,6.4,-1.5,522
+4,5,Somalia,709,-28.4,319.1,,,-27.7,-42.4,50.3,6.4,-1.5,547
+4,6.7,Northwest_Africa,714,-87.5,358.2,,,,,,,,7.5
+3,5.2,Northwest_Africa,714,-77.8,326.2,,,,,,,,8
+4,4.1,Northwest_Africa,714,-81.9,294.4,,,,,,,,13
+3,2.3,Northwest_Africa,714,-86.8,202.9,,,,,,,,13
+5,8,Northwest_Africa,714,-72,71.2,,,,,,,,81
+6,6.3,Northwest_Africa,714,-65.2,20.3,,,-65.3,25.8,33.6,26,2.3,152.5
+3,19.2,Northwest_Africa,714,-62.5,61.6,,,-61.3,65.8,33.6,26,2.3,160
+5,7.4,Northwest_Africa,714,-68.5,62.4,,,-67.3,67.4,33.6,26,2.3,185.5
+4,4.1,Northwest_Africa,714,-69.4,52,,,-68.5,57.7,33.6,26,2.3,187
+4,6.1,Northwest_Africa,714,-71.4,60.2,,,-70.2,66,33.6,26,2.3,187
+4,6.2,Northwest_Africa,714,-82.9,32.7,,,-82.4,48.7,33.6,26,2.3,193
+5,4.1,Northwest_Africa,714,-73,64.7,,,-71.7,70.6,33.6,26,2.3,200
+6,2.6,Northwest_Africa,714,-70.9,55.1,-76.2,78.9,-74.6,84.4,33.6,26,2.3,206.5
+5,6,Northwest_Africa,714,-29,60,-26.8,56.5,-25.7,58.6,33.6,26,2.3,273
+5,3.6,Northwest_Africa,714,-29.1,57.8,-26.3,54.2,-25.4,56.3,33.6,26,2.3,275
+5,2.8,Northwest_Africa,714,-38.5,57.5,-33.7,52.4,-32.8,54.8,33.6,26,2.3,286.5
+4,4.1,Northwest_Africa,714,-33.8,61.4,-29,55.5,-28,57.7,33.6,26,2.3,290
+5,3.5,Northwest_Africa,714,-28.7,55.8,-20.9,48.3,-20.1,50.3,33.6,26,2.3,307
+6,4.6,Northwest_Africa,714,-28.3,58.9,-21.1,51.5,-20.2,53.4,33.6,26,2.3,309
+5,2.6,Northwest_Africa,714,-32.8,55.7,-27.8,49.9,-27,52.1,33.6,26,2.3,310
+4,4.5,Northwest_Africa,714,-28.2,55.5,-20.3,47.8,-19.6,49.8,33.6,26,2.3,317
+7,5.3,Northwest_Africa,714,-26.6,44.7,-17.5,36.2,-17.1,38.1,33.6,26,2.3,320
+5,5.9,Northwest_Africa,714,-18.8,31.2,,,-18.6,33.2,33.6,26,2.3,348
+5,3.7,Northwest_Africa,714,-21,19.9,,,-21.2,21.9,33.6,26,2.3,365
+6,3.7,Northwest_Africa,714,-19.2,19.8,,,-19.4,21.8,33.6,26,2.3,365
+5,6.6,Northwest_Africa,714,-43.4,8.6,,,-43.9,11.7,33.6,26,2.3,409
+5,2.5,Northeast_Africa,715,-87.6,346.9,,,,,,,,1.5
+5,3.1,Northeast_Africa,715,-84.9,127.9,,,,,,,,3
+5,7.4,Northeast_Africa,715,-86.1,336.5,,,,,,,,4
+5,11.2,Northeast_Africa,715,-69,4,,,,,,,,11.5
+5,8.3,Northeast_Africa,715,-78.4,16.1,,,,,,,,11.5
+5,2.6,Northeast_Africa,715,-73,0.5,-83.2,341.7,,,,,,13.5
+4,12.7,Northeast_Africa,715,-83,13.3,,,,,,,,26.5
+6,6,Northeast_Africa,715,-79.6,332.2,-79.5,258.6,,,,,,29
+5,5.4,Northeast_Africa,715,-77.9,32.8,,,,,,,,30
+5,6.4,Northeast_Africa,715,-83.5,318.6,,,,,,,,37
+6,4,Northeast_Africa,715,-71,340,-76.5,308.4,,,,,,37.5
+6,4.2,Northeast_Africa,715,-78.1,342.8,,,,,,,,42.5
+6,6,Northeast_Africa,715,-68,338,-73,314.8,,,,,,44.5
+3,5.8,Northeast_Africa,715,-69.4,9.4,,,,,,,,44.5
+5,8.5,Northeast_Africa,715,-69.3,78.1,,,-69.4,78.3,39.9,-61.4,-0.2,93
+3,18.1,Northeast_Africa,715,-75.7,48.3,,,-75.9,48.4,40.1,-61.4,-0.2,94.5
+3,11.5,Northeast_Africa,715,-54.9,43.3,-59.7,27,-60.2,26.5,40.5,-61.4,-0.7,221.5
+4,5.5,Northeast_Africa,715,-54.5,45.8,,,-55,45.6,40.5,-61.4,-0.7,231
+5,3.8,Northeast_Africa,715,-59.3,34.1,,,-59.8,33.7,40.5,-61.4,-0.7,231
+4,6,Northeast_Africa,715,-40.8,71.3,,,-41.2,71.2,40.4,-61.4,-0.7,280
+3,10.8,Northeast_Africa,715,25.9,11.6,,,25.4,11.2,40.4,-61.4,-0.7,377
+6,9.2,Northeast_Africa,715,39.6,329.5,,,39.3,-30.6,40.4,-61.4,-0.7,463
+4,5.4,East_Gondwana,501,-39.2,105.6,,,-73.4,66.6,18.8,22.6,-39.2,64
+6,5.7,East_Gondwana,501,-39,100.8,,,-72,54.5,19,21.9,-40.2,65
+4,6.7,East_Gondwana,501,-40,96,,,-70.6,42.6,19.2,21.5,-40.7,65.5
+6,5.9,East_Gondwana,501,-38.4,102.4,,,-72.7,57.9,19.2,21.5,-40.7,65.5
+3,3.8,East_Gondwana,501,-34.5,103.6,,,-69.8,66.8,19.2,21.5,-40.7,65.5
+5,9.4,East_Gondwana,501,-37.2,100.5,,,-70.8,56.2,19.2,21.5,-40.7,65.5
+3,3.8,East_Gondwana,501,-39,99,,,-71.5,49.9,19.2,21.5,-40.7,65.5
+4,10.1,East_Gondwana,501,-34.6,94,,,-67.6,42.9,20.2,19.3,-43.8,69
+4,12,East_Gondwana,501,-21.6,119.4,,,-74.9,73.7,19.8,27.2,-59.2,88
+5,7.5,East_Gondwana,501,-14.2,117.8,,,-66.7,79.4,20.2,27.6,-58.5,91.2
+3,4,East_Gondwana,501,-3,118,,,-49.8,86.8,23.4,31.3,-53.8,116
+6,3.5,East_Gondwana,501,-7,117,,,-53.3,83.2,23.4,31.3,-53.8,116
+6,8.3,East_Gondwana,501,-9.3,124.8,,,-57.6,95.2,23.4,31.3,-53.8,116
+4,5.5,East_Gondwana,501,-6.5,120.2,,,-53.8,88.5,23.4,31.3,-53.8,116
+3,7,East_Gondwana,501,-16,121,,,-63,84.4,23.5,31.4,-53.7,116.5
+5,2.4,East_Gondwana,501,-9.4,116.6,,,-55.1,81.1,23.5,31.5,-53.6,117
+5,4.6,East_Gondwana,501,-10.1,130.1,-2.4,118.5,-45.1,71.3,29.8,42.1,-60.5,206.5
+4,4.6,East_Gondwana,501,7.3,124.3,12.2,110.9,-28.7,73.5,29.8,42.1,-60.5,243
+5,6,East_Gondwana,501,7.5,120.5,13.4,109.1,-26.7,72.5,29.8,42.1,-60.5,248.5
+6,4.3,East_Gondwana,501,2.2,125.8,,,-44.6,83.5,29.8,42.1,-60.5,250.5
+3,1.8,East_Gondwana,501,4.1,102.8,10.2,94.5,-20.3,58.1,29.8,42.1,-60.5,250.5
+6,6.5,East_Gondwana,501,4,129,9.5,115.1,-33.2,75.8,29.8,42.1,-60.5,250.5
+3,12.1,East_Gondwana,501,18.1,111,,,-23.8,77,29.8,42.1,-60.5,289.5
+5,4.1,Arabia,503,-82.4,62.2,,,-82.9,57.8,36.5,18,-0.8,2.5
+5,3.4,Arabia,503,-76,13.5,,,-75.2,-3.1,36.5,18,-4.5,19
+5,3.6,Arabia,503,-74.2,69.1,,,-77.7,46.7,35.7,20.4,-7.1,29
+3,7.2,Arabia,503,-25.6,64,-13.4,58.4,-17.7,51.5,37.1,17.1,-8.9,306.5
+5,4.4,Madagascar,702,-74,43.7,,,,,,,,86.5
+4,7.6,Madagascar,702,-64,63,,,,,,,,87
+4,4.9,Madagascar,702,-66.1,49.7,,,,,,,,87
+4,4.4,Madagascar,702,-65.8,35.6,,,,,,,,87
+4,8.9,Madagascar,702,-73.7,73.1,,,,,,,,87
+5,4.3,Madagascar,702,-65.5,38,,,,,,,,87
+4,6.9,Madagascar,702,-70.3,63.1,,,,,,,,87
+5,2.4,Madagascar,702,-76.8,68.2,,,,,,,,87
+5,10.7,Madagascar,702,-66.7,43.5,,,,,,,,91
+3,5.9,Madagascar,702,-74,97.1,-65.2,70.8,-50.9,59.5,14.8,137.5,-15.4,206.5
+3,7.6,Madagascar,702,-76.7,110.8,-68.4,73.5,-54.3,60.4,14.8,137.5,-15.4,250.5
+3,9.5,Madagascar,702,-51.3,72.6,-42.5,60.2,-27.8,54.4,14.8,137.5,-15.4,305
+5,11,Madagascar,702,-6.8,1,,,3.1,-2.6,14.8,137.5,-15.4,509
+4,14,Madagascar,702,-6.8,352.7,,,1.6,-10.7,14.8,137.5,-15.4,521
+5,1.8,Australia,801,-78.4,103.1,-70.3,148.6,-80.2,155.3,12,48.8,-10.4,17.5
+5,3.6,Australia,801,-70.5,125.6,,,-80.2,19,14.1,57,-24.4,53
+5,1.8,Australia,801,-48.9,148.7,-45.1,146.8,-56.4,93.9,17.2,103.2,-36,112
+5,3.7,Australia,801,-53,179.6,,,-59.2,64.5,19.5,117.8,-56.2,168
+4,6,Australia,801,-46.1,175.2,,,-57.3,78.3,19.5,117.8,-56.2,180
+5,2.9,Australia,801,-50.7,174.5,,,-56.6,69.9,19.5,117.8,-56.2,183
+6,5.1,Australia,801,-63.8,124.5,,,-31.5,48.8,19.5,117.8,-56.2,321
+6,6,Australia,801,-47.1,41,-49.3,54.7,2.2,28.5,19.5,117.8,-56.2,367
+6,7.8,Australia,801,-49.1,38,,,3.6,17.7,19.5,117.8,-56.2,370
+6,15.2,Australia,801,-62,23.2,,,-10,10.7,19.5,117.8,-56.2,370
+5,5,Australia,801,-63,38.6,,,-10.3,18,19.5,117.8,-56.2,376
+5,3,Australia,801,-26.7,33.7,-30.4,46,21.9,25.2,19.5,117.8,-56.2,465
+4,13,Australia,801,-13,25,,,38.3,1.5,19.5,117.8,-56.2,493
+6,7.4,Australia,801,3.1,54.1,,,52.7,45,19.5,117.8,-56.2,495
+3,3.8,Australia,801,-37.5,34.4,,,15.2,14.8,19.5,117.8,-56.2,500
+6,10,Australia,801,-19.3,39.1,,,33.4,19,19.5,117.8,-56.2,510
+3,10.1,Australia,801,-31.4,26.9,-31.3,26.7,20.6,7.4,19.5,117.8,-56.2,510
+3,10.4,Australia,801,-38.3,24.5,,,13.4,6.9,19.5,117.8,-56.2,510
+3,12.3,Australia,801,-33.8,15.1,-32.8,13.9,16.4,-3.2,19.5,117.8,-56.2,515
+7,6.7,Australia,801,-43.2,339.9,,,-6.2,-20.8,19.5,117.8,-56.2,522
+6,14.4,Australia,801,-37.4,20.1,,,13.5,3.2,19.5,117.8,-56.2,523
+6,7.3,Australia,801,-32.7,11.5,-28.7,5.7,17.4,-11.7,19.5,117.8,-56.2,525
+6,4.1,Australia,801,-46.6,337.3,-37.4,331.7,-6.5,-29.4,19.5,117.8,-56.2,534
+5,11.4,Australia,801,-21.3,14.9,,,27.5,-6.6,19.5,117.8,-56.2,535
+6,16,Australia,801,-33,328,-20.2,326.5,1.6,-45.3,19.5,117.8,-56.2,550
+5,4,East_Antarctica,802,-85.5,143.6,,,-85.5,141.7,8.3,-49.4,0.2,1
+5,6.3,East_Antarctica,802,-87.3,137.3,,,-87.3,130.8,8.2,-49.4,0.3,2
+5,7,East_Antarctica,802,-86.4,28.5,,,-85.7,30.7,8.2,-49.4,0.8,5
+5,2.3,East_Antarctica,802,-85.5,9.3,,,-81.7,25.3,11.6,-48.2,4.2,27
+5,4.4,East_Antarctica,802,-51.4,203.4,,,-55,98.6,10.5,148.8,-58.2,164
+4,3.4,East_Antarctica,802,-47.8,225.5,,,-69.3,91.3,10.5,148.8,-58.2,183
+3,3.3,East_Antarctica,802,-57.8,224.3,,,-61.5,75.7,10.5,148.8,-58.2,183
+5,2.4,East_Antarctica,802,-45.3,208,,,-59.4,108.3,10.5,148.8,-58.2,183
+5,10.2,East_Antarctica,802,-50.5,211.4,,,-60,97.1,10.5,148.8,-58.2,183
+4,6.9,East_Antarctica,802,-44.1,231.5,,,-74.9,91.3,10.5,148.8,-58.2,193
+5,3.8,East_Antarctica,802,-41.8,226.5,,,-73.2,105.8,10.5,148.8,-58.2,195
+5,5.2,East_Antarctica,802,-2.5,23.8,,,37.8,-2.5,10.5,148.8,-58.2,471.5
+5,7.6,East_Antarctica,802,-3.5,22.7,,,36.3,-2.7,10.5,148.8,-58.2,475
+5,7.2,East_Antarctica,802,-11,21,,,29.6,2,10.5,148.8,-58.2,479
+4,10.9,East_Antarctica,802,-9.3,26.7,,,34.5,5.7,10.5,148.8,-58.2,484
+4,12,East_Antarctica,802,-7,21.4,,,32.8,-0.9,10.5,148.8,-58.2,499
+4,8.1,East_Antarctica,802,-5.4,18.5,,,32,-4.6,10.5,148.8,-58.2,500
+3,4.5,East_Antarctica,802,-28.5,9.5,,,9.5,6.1,10.5,148.8,-58.2,515
diff --git a/examples/Torsvik_2012/data.jld2 b/examples/Torsvik_2012/data.jld2
new file mode 100644
index 0000000..eeccae0
Binary files /dev/null and b/examples/Torsvik_2012/data.jld2 differ
diff --git a/examples/Torsvik_2012/plots/plot_results.jl b/examples/Torsvik_2012/plots/plot_results.jl
new file mode 100644
index 0000000..eaa314a
--- /dev/null
+++ b/examples/Torsvik_2012/plots/plot_results.jl
@@ -0,0 +1,140 @@
+using Pkg; Pkg.activate(".")
+
+using SphereUDE
+using Serialization
+using Plots 
+using Plots.PlotMeasures
+using JLD2
+using LinearAlgebra
+
+JLD2.@load "examples/Torsvik_2012/results/results_dict.jld2" 
+
+data_directions_sph = cart2sph(results_dict["directions"], radians=false)
+fit_directions_sph = cart2sph(results_dict["fit_directions"], radians=false)
+
+data_lat = data_directions_sph[1,:]
+fit_lat  = fit_directions_sph[1,:]
+
+data_lon = data_directions_sph[2,:]
+fit_lon  = fit_directions_sph[2,:]
+
+α_error = 140 ./ sqrt.(results_dict["kappas"])
+
+# Latitude Plot 
+
+Plots.scatter(results_dict["times"], data_lat, label="Paleopole Latitudes", yerr=α_error, c=:lightsteelblue2, markersize=5)
+plot_latitude = Plots.plot!(results_dict["fit_times"], fit_lat, label="Estimated APWP using SphereUDE", 
+                    xlabel="Age (Myr)", yticks=[-90, -60, -30, 0, 30, 60],
+                    ylabel="Latitude (degrees)", ylims=(-90,60), lw = 4, c=:brown,
+                    legend=:topleft)
+plot!(fontfamily="Computer Modern",
+     #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    legend=true,
+    margin= 7mm,
+    size=(1200,500),
+    dpi=600)
+
+Plots.savefig(plot_latitude, "examples/Torsvik_2012/plots/latitude.pdf")
+
+### Longitude Plot 
+
+α_error_lon = α_error ./ cos.(π ./ 180. .* data_lat)
+
+Plots.scatter(results_dict["times"], data_lon, label="Paleopole Longitudes", yerr=α_error_lon, c=:lightsteelblue2, markersize=5)
+plot_longitude = Plots.plot!(results_dict["fit_times"], fit_lon, label="Estimated APWP using SphereUDE", 
+                    xlabel="Age (Myr)", yticks=[-180, -120, -60 , 0, 60, 120, 180],
+                    ylabel="Longitude (degrees)", ylims=(-180,180), lw = 4, c=:brown,
+                    legend=:bottomright)
+plot!(fontfamily="Computer Modern",
+     #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    margin= 7mm,
+    size=(1200,500),
+    dpi=600)
+
+Plots.savefig(plot_longitude, "examples/Torsvik_2012/plots/longitude.pdf")
+
+### Lat and long combined
+
+combo_plot = plot(plot_latitude, plot_longitude, layout = (2, 1))
+plot!(fontfamily="Computer Modern",
+     #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    margin= 7mm,
+    size=(1200,800),
+    dpi=600)
+Plots.savefig(combo_plot, "examples/Torsvik_2012/plots/latitude_longitude.pdf")
+
+
+### Angular velocity Plot
+
+angular_velocity = mapslices(x -> norm(x), results_dict["fit_rotations"], dims=1)[:]
+angular_velocity_path = [norm(cross(results_dict["fit_directions"][:,i], results_dict["fit_rotations"][:,i] )) for i in axes(results_dict["fit_rotations"],2)]
+
+plot_angular = Plots.plot(results_dict["fit_times"], angular_velocity, label="Maximum total angular velocity", 
+                    xlabel="Age (Myr)", 
+                    ylabel="Angular velocity (degrees/My)", lw = 5, c=:darkgreen,
+                    legend=:topleft)
+plot!(results_dict["fit_times"], angular_velocity_path, label="Pole angular velocity", lw = 4, c=:lightseagreen, ls=:dot)
+plot!(fontfamily="Computer Modern",
+    #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    margin= 10mm,
+    size=(1200,500),
+    dpi=300)
+
+
+Plots.savefig(plot_angular, "examples/Torsvik_2012/plots/angular.pdf")
+
+
+### Loss function
+
+losses = results_dict["losses"]
+
+plot_loss = Plots.plot(1:length(losses), losses, label="Loss Function", 
+                    xlabel="Epoch", 
+                    ylabel="Loss", lw = 5, c=:indigo,
+                    yaxis=:log,
+                    # yticks=[1,10,100],
+                    xlimits=(0,10000),
+                    ylimits=(10, 1000),
+                    xticks=[0,2500,5000,7500,10000],
+                    legend=:topright)
+
+vspan!(plot_loss, [0,5000], color = :navajowhite3, alpha = 0.2, labels = "ADAM");
+vspan!(plot_loss, [5000,10000], color = :navajowhite4, alpha = 0.2, labels = "BFGS");
+
+plot!(fontfamily="Computer Modern",
+    #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    margin= 10mm,
+    size=(700,500),
+    dpi=300)
+
+Plots.savefig(plot_loss, "examples/Torsvik_2012/plots/loss.pdf")
diff --git a/examples/Torsvik_2012/results.jld2 b/examples/Torsvik_2012/results.jld2
new file mode 100644
index 0000000..40676b8
Binary files /dev/null and b/examples/Torsvik_2012/results.jld2 differ
diff --git a/examples/Torsvik_2012/results/_results_dict.jld2 b/examples/Torsvik_2012/results/_results_dict.jld2
new file mode 100644
index 0000000..991dfbd
Binary files /dev/null and b/examples/Torsvik_2012/results/_results_dict.jld2 differ
diff --git a/examples/Torsvik_2012/results/data.jld2 b/examples/Torsvik_2012/results/data.jld2
new file mode 100644
index 0000000..eeccae0
Binary files /dev/null and b/examples/Torsvik_2012/results/data.jld2 differ
diff --git a/examples/Torsvik_2012/results/results.jld2 b/examples/Torsvik_2012/results/results.jld2
new file mode 100644
index 0000000..5e1e7e5
Binary files /dev/null and b/examples/Torsvik_2012/results/results.jld2 differ
diff --git a/examples/Torsvik_2012/results/results_dict.jld2 b/examples/Torsvik_2012/results/results_dict.jld2
new file mode 100644
index 0000000..eb9d99b
Binary files /dev/null and b/examples/Torsvik_2012/results/results_dict.jld2 differ
diff --git a/examples/double_rotation/2rotations.ipynb b/examples/double_rotation/2rotations.ipynb
index a3fefe3..2ea732d 100644
--- a/examples/double_rotation/2rotations.ipynb
+++ b/examples/double_rotation/2rotations.ipynb
@@ -29,7 +29,7 @@
     "using SciMLSensitivity\n",
     "using Optimization, OptimizationOptimisers, OptimizationOptimJL\n",
     "\n",
-    "using SphereFit"
+    "using SphereUDE"
    ]
   },
   {
diff --git a/examples/double_rotation/plot_L.pdf b/examples/double_rotation/plot_L.pdf
deleted file mode 100644
index 1b215b3..0000000
Binary files a/examples/double_rotation/plot_L.pdf and /dev/null differ
diff --git a/examples/double_rotation/plot_sphere.pdf b/examples/double_rotation/plot_sphere.pdf
deleted file mode 100644
index a199811..0000000
Binary files a/examples/double_rotation/plot_sphere.pdf and /dev/null differ
diff --git a/examples/multi_rotation/plot_sphere.pdf b/examples/multi_rotation/plot_sphere.pdf
index e59ff07..83d32f7 100644
Binary files a/examples/multi_rotation/plot_sphere.pdf and b/examples/multi_rotation/plot_sphere.pdf differ
diff --git a/examples/reg_gradient.md b/examples/reg_gradient.md
new file mode 100644
index 0000000..66fd43f
--- /dev/null
+++ b/examples/reg_gradient.md
@@ -0,0 +1,177 @@
+How to handle more than one AD method in SciML applications
+
+Hi all, 
+
+I am implementing regularization penalties inside Universal Differential Equations (also applicable to Physics-Informed neural networks) where I need to differentiate a (loss) function that includes in its calculation a gradient. For an UDE, consider the case where the outputs of the neural network are differentiated with respect to the input layer and added to the loss function (this will represent a physical derivative that I am interested in constraining for some application). A second (computational) derivative with respect to the weights of the NN needs to be computed on top of this for optimizing the weights of the NN. 
+
+Consider the simple example based on the [Lotka-Volterra equations from the SciML tutorial](https://docs.sciml.ai/Overview/dev/showcase/missing_physics/) with the following (slightly modified) loss function (MWE with all code is provided at the end of this post):
+```julia 
+function loss(θ)
+    # Empirical error
+    X̂ = predict(θ)
+    l_emp = mean(abs2, Xₙ .- X̂)
+    # Regularization based on first derivatives 
+    steps_reg = collect(0.0:0.1:10.0)
+    dUdx = map(x -> Zygote.jacobian(first ∘ U, [x, 1.0], θ, st)[1], steps_reg)
+    norm_dUdx = norm.(dUdx).^2.0
+    l_reg = sum(norm_dUdx) 
+    return l_emp + l_reg
+end
+```
+The loss function here is the combination of the empirical error and a regularization term involving derivatives of the neural network. This can be easily evaluated, so the following works
+```julia 
+loss(p)
+# Return: 116666.31176556791
+```
+However, when applying the `solve` step to optimize the weights of the neural network I obtain `ERROR: LoadError: Mutating arrays is not supported`, but I am somehow confident the problem is coming from mutating arrays in my code. For example. if I change the loss function to not compute this extra gradient, everything works fine:
+```julia
+function loss(θ)
+    # Empirical error
+    X̂ = predict(θ)
+    l_emp = mean(abs2, Xₙ .- X̂)
+    # regularization 
+    steps_reg = collect(0.0:0.1:10.0)
+    # dUdx = map(x -> Zygote.jacobian(first ∘ U, [x, 1.0], θ, st)[1], steps_reg) # previous code
+    dUdx = map(x -> U([x, 1.0], θ, st)[1], steps_reg) # new code
+    norm_dUdx = norm.(dUdx).^2.0
+    l_reg = sum(norm_dUdx) 
+    return l_emp + l_reg
+end
+```
+A second implementation of this uses finite differences to compute the derivative involved in the regularization. 
+```julia 
+function loss(θ)
+    # Empirical error
+    X̂ = predict(θ)
+    l_emp = mean(abs2, Xₙ .- X̂)
+    # Regularization computed with finite differences  
+    steps_reg = collect(0.0:0.1:10.0)
+    Ux = map(x -> (first ∘ U)([x, 1.0], θ, st), steps_reg)
+    dUdx = diff(Ux) ./ diff(steps_reg)
+    l_reg = sum(norm.(dUdx).^2.0)    
+    return l_emp + l_reg
+end
+```
+This of course works, although I noticed that it is slower than expected (especially compared to the adjoint differentiation involved in the numerical solver). However, this is far from being the ideal solution to the problem, since I would like to be able to compute derivatives inside my loss function with some AD tool.
+
+*The big picture problem.* I think the really interesting thing here is how different differentiation methods interact with each other. This is also relevant for physics-informed neural networks, where most implementations I have seen rely purely on reverse AD to compute higher-order derivatives (both with respect to input layer and weights). However, this is quite inefficient and scales exponentially with the order of the derivative. On the other side, I noticed that `NeuralPDE.jl` instead uses the methods of lines to discretize derivatives in space (code [here](https://github.com/SciML/NeuralPDE.jl/blob/ce31a5ee84fd5ecd6c1532b4fc175d4b3eddd11a/src/pinn_types.jl#L379-L414)), which is quite the same I am doing with this example. 
+
+*The ideal solution.* I would like to have something that allows me to compute (physical) derivatives using some form of AD (for derivatives with respect to input layers of a NN, probably a forward method) and then be able to perform gradient-based optimization on top of it using reverse AD. This will push very interesting SciML applications with physical regularization (I have a few examples in my own research) and it is becoming increasingly important for deep learning methods used to solve differential equations. 
+
+I would appreciate if someone had some insights on this, both in 
+1) Do you have an idea how to compute the gradients in the loss function using some form of AD methods for the regularization; and
+2) Big picture perspectives of the integration between different AD tools in real SciML problems. 
+Thanks! 
+
+Complete minimal working example:
+```julia
+using Pkg; Pkg.activate(".")
+
+# SciML Tools
+using OrdinaryDiffEq, SciMLSensitivity
+using Optimization, OptimizationOptimisers, OptimizationOptimJL
+using LinearAlgebra, Statistics
+using ComponentArrays, Lux, Zygote #, StableRNGs
+
+# Set a random seed for reproducible behaviour
+using Random
+rng = Random.default_rng()
+Random.seed!(rng, 666)
+
+function lotka!(du, u, p, t)
+    α, β, γ, δ = p
+    du[1] = α * u[1] - β * u[2] * u[1]
+    du[2] = γ * u[1] * u[2] - δ * u[2]
+end
+
+# Define the experimental parameter
+tspan = (0.0, 5.0)
+u0 = 5.0f0 * rand(rng, 2)
+p_ = [1.3, 0.9, 0.8, 1.8]
+prob = ODEProblem(lotka!, u0, tspan, p_)
+solution = solve(prob, Vern7(), abstol = 1e-12, reltol = 1e-12, saveat = 0.25)
+
+# Add noise in terms of the mean
+X = Array(solution)
+t = solution.t
+
+x̄ = mean(X, dims = 2)
+noise_magnitude = 5e-3
+Xₙ = X .+ (noise_magnitude * x̄) .* randn(rng, eltype(X), size(X))
+
+rbf(x) = exp.(-(x .^ 2))
+# Multilayer FeedForward
+const U = Lux.Chain(Lux.Dense(2, 5, rbf), Lux.Dense(5, 5, rbf), Lux.Dense(5, 5, rbf),
+              Lux.Dense(5, 2))
+# Get the initial parameters and state variables of the model
+p, st = Lux.setup(rng, U)
+const _st = st
+
+# Define the hybrid model
+function ude_dynamics!(du, u, p, t, p_true)
+    û = U(u, p, _st)[1] # Network prediction
+    du[1] = p_true[1] * u[1] + û[1]
+    du[2] = -p_true[4] * u[2] + û[2]
+end
+
+# Closure with the known parameter
+nn_dynamics!(du, u, p, t) = ude_dynamics!(du, u, p, t, p_)
+# Define the problem
+prob_nn = ODEProblem(nn_dynamics!, Xₙ[:, 1], tspan, p)
+
+function predict(θ, X = Xₙ[:, 1], T = t)
+    _prob = remake(prob_nn, u0 = X, tspan = (T[1], T[end]), p = θ)
+    Array(solve(_prob, Vern7(), saveat = T,
+                abstol = 1e-6, reltol = 1e-6,
+                sensealg=QuadratureAdjoint(autojacvec=ReverseDiffVJP(true))))
+end
+
+"""
+Works but it is clearly a little bit slow! 
+"""
+# function loss(θ)
+#     # Empirical error
+#     X̂ = predict(θ)
+#     l_emp = mean(abs2, Xₙ .- X̂)
+#     # regularization 
+#     steps_reg = collect(0.0:0.1:10.0)
+#     Ux = map(x -> (first ∘ U)([x, 1.0], θ, st), steps_reg)
+#     dUdx = diff(Ux) ./ diff(steps_reg)
+#     l_reg = sum(norm.(dUdx).^2.0)    
+#     return l_emp + l_reg
+# end
+
+"""
+Need to solve problem with mutating arrays!
+"""
+function loss(θ)
+    # Empirical error
+    X̂ = predict(θ)
+    l_emp = mean(abs2, Xₙ .- X̂)
+    # regularization 
+    steps_reg = collect(0.0:0.1:10.0)
+    dUdx = map(x -> Zygote.jacobian(first ∘ U, [x, 1.0], θ, st)[1], steps_reg)
+    # dUdx = map(x -> U([x, 1.0], θ, st)[1], steps_reg)
+    norm_dUdx = norm.(dUdx).^2.0
+    l_reg = sum(norm_dUdx) 
+    return l_emp + l_reg
+end
+
+losses = Float64[]
+
+callback = function (p, l)
+    push!(losses, l)
+    if length(losses) % 10 == 0
+        println("Current loss after $(length(losses)) iterations: $(losses[end])")
+    end
+    return false
+end
+
+adtype = Optimization.AutoZygote()
+# adtype = Optimization.AutoReverseDiff()
+optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
+optprob = Optimization.OptimizationProblem(optf, ComponentVector{Float64}(p))
+
+res1 = Optimization.solve(optprob, ADAM(), callback = callback, maxiters = 500)
+println("Training loss after $(length(losses)) iterations: $(losses[end])")
+```
\ No newline at end of file
diff --git a/examples/single_rotation_time/simulate.jl b/examples/single_rotation_time/simulate.jl
new file mode 100644
index 0000000..f8544ff
--- /dev/null
+++ b/examples/single_rotation_time/simulate.jl
@@ -0,0 +1,36 @@
+using Pkg; Pkg.activate(".")
+
+using SphereUDE
+
+# Random seed
+using Random, Distributions
+rng = Random.default_rng()
+Random.seed!(rng, 616)
+
+# Total time simulation
+tspan = [0, 160.0]
+# Number of sample points
+N_samples = 100
+
+# The original time series consist in an uniform grid on time
+times_true = collect(range(tspan[1], tspan[2], N_samples+1))[2:end]
+# Amoung of noise in timeseroes 
+σ = 0.1
+times_samples = times_true .+ rand(sampler(Normal(0, σ)), N_samples)
+
+
+
+
+# Expected maximum angular deviation in one unit of time (degrees)
+Δω₀ = 1.0   
+# Angular velocity 
+ω₀ = Δω₀ * π / 180.0
+# Change point
+τ₀ = 65.0
+# Angular momentum
+L0 = ω₀    .* [1.0, 0.0, 0.0]
+L1 = 0.5ω₀ .* [0.0, 1/sqrt(2), 1/sqrt(2)]
+
+# Solver tolerances 
+reltol = 1e-7
+abstol = 1e-7
\ No newline at end of file
diff --git a/examples/test_CS.jl b/examples/test_CS.jl
new file mode 100644
index 0000000..78f0198
--- /dev/null
+++ b/examples/test_CS.jl
@@ -0,0 +1,98 @@
+using Pkg; Pkg.activate(".")
+using Revise 
+
+using LinearAlgebra, Statistics, Distributions 
+using OrdinaryDiffEq
+using SciMLSensitivity
+using Optimization, OptimizationOptimisers, OptimizationOptimJL
+using ComponentArrays, Lux, Zygote
+
+##############################################################
+###############  Simulation of Simple Example ################
+##############################################################
+
+# Random seed
+using Random
+rng = Random.default_rng()
+Random.seed!(rng, 666)
+
+function lotka!(du, u, p, t)
+    α, β, γ, δ = p
+    du[1] = α * u[1] - β * u[2] * u[1]
+    du[2] = γ * u[1] * u[2] - δ * u[2]
+end
+
+# Define the experimental parameter
+tspan = (0.0, 5.0)
+u0 = 5.0f0 * rand(rng, 2)
+p_ = [1.3, 0.9, 0.8, 1.8]
+prob = ODEProblem(lotka!, u0, tspan, p_)
+solution = solve(prob, Vern7(), abstol = 1e-12, reltol = 1e-12, saveat = 0.25)
+
+# Add noise in terms of the mean
+X = Array(solution)
+t = solution.t
+
+x̄ = mean(X, dims = 2)
+noise_magnitude = 5e-3
+Xₙ = X .+ (noise_magnitude * x̄) .* randn(rng, eltype(X), size(X))
+
+# plot(solution, alpha = 0.75, color = :black, label = ["True Data" nothing])
+# scatter!(t, transpose(Xₙ), color = :red, label = ["Noisy Data" nothing])
+
+rbf(x) = exp.(-(x .^ 2))
+
+# Multilayer FeedForward
+const U = Lux.Chain(Lux.Dense(2, 5, rbf), Lux.Dense(5, 5, rbf), Lux.Dense(5, 5, rbf),
+              Lux.Dense(5, 2))
+# Get the initial parameters and state variables of the model
+p, st = Lux.setup(rng, U)
+const _st = st
+
+# Define the hybrid model
+function ude_dynamics!(du, u, p, t, p_true)
+    û = U(u, p, _st)[1] # Network prediction
+    du[1] = p_true[1] * u[1] + û[1]
+    du[2] = -p_true[4] * u[2] + û[2]
+end
+
+# Closure with the known parameter
+nn_dynamics!(du, u, p, t) = ude_dynamics!(du, u, p, t, p_)
+# Define the problem
+prob_nn = ODEProblem(nn_dynamics!, Xₙ[:, 1], tspan, p)
+
+function predict(θ, T = t)
+    _prob = remake(prob_nn, u0 = θ.u0, tspan = (T[1], T[end]), p = θ.p)
+    Array(solve(_prob, Vern7(), saveat = T,
+                abstol = 1e-6, reltol = 1e-6,
+                sensealg=QuadratureAdjoint(autojacvec=ReverseDiffVJP(true))))
+end
+
+function loss(θ)
+    # Empirical error
+    X̂ = predict(θ)
+    l_emp = mean(abs2, Xₙ .- X̂)
+    # l_emp = imag(1 + mean(Xₙ .- X̂) * im)
+    return l_emp
+end
+
+losses = Float64[]
+
+callback = function (p, l)
+    push!(losses, l)
+    if length(losses) % 100 == 0
+        println("Current loss after $(length(losses)) iterations: $(losses[end])")
+    end
+    return false
+end
+
+β = ComponentVector(p=p, u0=u0)
+
+adtype = Optimization.AutoZygote()
+# adtype = Optimization.AutoReverseDiff()
+# adtype = Optimization.AutoForwardDiff()
+optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
+optprob = Optimization.OptimizationProblem(optf, β)
+
+res1 = Optimization.solve(optprob, ADAM(), callback = callback, maxiters = 500)
+println("Training loss after $(length(losses)) iterations: $(losses[end])")
diff --git a/examples/working-example.jl b/examples/working-example.jl
index d763d48..dd8ac89 100644
--- a/examples/working-example.jl
+++ b/examples/working-example.jl
@@ -20,6 +20,8 @@ Random.seed!(rng, 666)
 # using TimerOutputs
 # const to = TimerOutput()
 
+# Testing new solution
+# using BatchedRoutines
 
 function lotka!(du, u, p, t)
     α, β, γ, δ = p
@@ -47,9 +49,14 @@ Xₙ = X .+ (noise_magnitude * x̄) .* randn(rng, eltype(X), size(X))
 
 rbf(x) = exp.(-(x .^ 2))
 
+
 # Multilayer FeedForward
-const U = Lux.Chain(Lux.Dense(2, 5, rbf), Lux.Dense(5, 5, rbf), Lux.Dense(5, 5, rbf),
-              Lux.Dense(5, 2))
+const U = Lux.Chain(
+                Lux.Dense(2, 5, rbf), 
+                Lux.Dense(5, 5, rbf), 
+                Lux.Dense(5, 5, rbf),
+                Lux.Dense(5, 2)
+                )
 # Get the initial parameters and state variables of the model
 p, st = Lux.setup(rng, U)
 const _st = st
@@ -73,6 +80,41 @@ function predict(θ, X = Xₙ[:, 1], T = t)
                 sensealg=QuadratureAdjoint(autojacvec=ReverseDiffVJP(true))))
 end
 
+
+# Code comming from Lux example on nested AD
+
+function loss_function1(model, x, ps, st, y)
+    # Make it a stateful layer
+    smodel = StatefulLuxLayer{true}(model, ps, st)
+    ŷ = smodel(x)
+    loss_emp = sum(abs2, ŷ .- y)
+    # You can use `Zygote.jacobian` as well but ForwardDiff tends to be more efficient here
+    J = ForwardDiff.jacobian(smodel, x)
+    loss_reg = abs2(norm(J .* 0.01f0))
+    return loss_emp + loss_reg
+end
+
+# Using Batchnorm to show that it is possible
+model = Chain(Dense(2 => 4, tanh), BatchNorm(4), Dense(4 => 2))
+ps, st = Lux.setup(StableRNG(0), model)
+x = randn(StableRNG(0), Float32, 2, 10)
+y = randn(StableRNG(11), Float32, 2, 10)
+
+loss_function1(model, x, ps, st, y)
+
+
+
+
+
+
+
+
+
+
+
+
+# From here is the old code for test
+
 # function loss(θ)
 #     # Empirical error
 #     X̂ = predict(θ)
@@ -104,7 +146,13 @@ function loss(θ)
     l_emp = mean(abs2, Xₙ .- X̂)
     # regularization 
     steps_reg = collect(0.0:0.1:10.0)
-    dUdx = map(x -> Zygote.jacobian(first ∘ U, [x, 1.0], θ, st)[1], steps_reg)
+    # dUdx = map(x -> Zygote.jacobian(first ∘ U, [x, 1.0], θ, st)[1], steps_reg)
+    # dUdx = map(x -> ForwardDiff.jacobian(x -> U([x[1], 1.0], θ, st)[1], [x]), steps_reg)
+    # @show batched_jacobian(AutoForwardDiff(), x -> U([x[1], 1.0], θ, st)[1], [1.0])
+    dUdx = map(x -> batched_jacobian(AutoForwardDiff(), x -> U([x[1], 1.0], θ, st)[1], [1.0])[1], steps_reg) # this works with adtype=AutoForwardDiff
+    # dUdx = map(x -> batched_jacobian(AutoReverseDiff(), x -> U([x[1], 1.0], θ, st)[1], [1.0])[1], steps_reg)
+    # dUdx = map(x -> batched_jacobian(AutoForwardDiff(), ), steps_reg)
+    # @show dUdx
     norm_dUdx = norm.(dUdx).^2.0
     l_reg = sum(norm_dUdx) 
     return l_emp + l_reg
@@ -120,7 +168,9 @@ callback = function (p, l)
     return false
 end
 
-adtype = Optimization.AutoZygote()
+# adtype = Optimization.AutoZygote()
+# adtype = Optimization.AutoReverseDiff()
+adtype = Optimization.AutoForwardDiff()
 optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
 optprob = Optimization.OptimizationProblem(optf, ComponentVector{Float64}(p))
 
diff --git a/file.dat b/file.dat
new file mode 100644
index 0000000..6fee52f
Binary files /dev/null and b/file.dat differ
diff --git a/magnus_test.jl b/magnus_test.jl
new file mode 100644
index 0000000..25df2e9
--- /dev/null
+++ b/magnus_test.jl
@@ -0,0 +1,23 @@
+using OrdinaryDiffEqCore, OrdinaryDiffEqTsit5
+using OrdinaryDiffEqLinear
+using BenchmarkTools
+
+function update_func(A, u, p, t)
+    A[1, 1] = cos(t)
+    A[2, 1] = sin(t)
+    A[1, 2] = -sin(t)
+    A[2, 2] = cos(t)
+end
+
+A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+prob = ODEProblem(A, ones(2), (1.0, 6.0))
+
+@btime solve(prob, Tsit5(), dt=0.1, adaptive=false);
+# 38.446 μs (487 allocations: 41.09 KiB)
+
+@btime solve(prob, MagnusGL6(), dt=0.1);
+# 208.365 μs (3481 allocations: 335.88 KiB)
+@btime solve(prob, MagnusGL4(), dt=0.1);
+# 93.735 μs (981 allocations: 94.47 KiB)
+@btime solve(prob, MagnusLeapfrog(), dt=0.1);
+# 67.176 μs (532 allocations: 45.33 KiB)
\ No newline at end of file
diff --git a/src/train_test.jl b/src/train_test.jl
new file mode 100644
index 0000000..f23494b
--- /dev/null
+++ b/src/train_test.jl
@@ -0,0 +1,133 @@
+export train
+
+function get_NN(params, rng, θ_trained)
+    # Define neural network 
+    U = Lux.Chain(
+        Lux.Dense(1,  5,  relu_cap), # explore discontinuity function for activation
+        Lux.Dense(5,  10, relu_cap), 
+        Lux.Dense(10, 5,  relu_cap), 
+        Lux.Dense(5,  3,  x->sigmoid_cap(x; ω₀=params.ωmax))
+    )
+    θ, st = Lux.setup(rng, U)
+    return U, θ, st
+end
+
+function train(data::AD,
+               params::AP,
+               rng, 
+               θ_trained=[]) where{AD <: AbstractData, AP <: AbstractParameters}
+
+    U, θ, st = get_NN(params, rng, θ_trained)
+
+    # one option is to restrict where the NN is evaluated to discrete t to 
+    # generate piece-wise dynamics.
+
+    function ude_rotation!(du, u, p, t)
+        # Angular momentum given by network prediction
+        L = U([t], p, st)[1]
+        du .= cross(L, u)
+    end
+
+    prob_nn = ODEProblem(ude_rotation!, params.u0, [params.tmin, params.tmax], θ)
+
+    function predict(θ::ComponentVector; u0=params.u0, T=data.times) 
+        _prob = remake(prob_nn, u0=u0, 
+                       tspan=(min(T[1], params.tmin), max(T[end], params.tmax)), 
+                       p = θ)
+        Array(solve(_prob, params.solver, saveat=T,
+                    abstol=params.abstol, reltol=params.reltol,
+                    sensealg=params.sensealg))
+    end
+
+    function loss(θ::ComponentVector)
+        u_ = predict(θ)
+        # Empirical error
+        # l_emp = mean(abs2, u_ .- data.directions)
+        if isnothing(data.kappas)
+            # The 3 is needed since the mean is computen on a 3xN matrix
+            l_emp = 1 - 3.0 * mean(u_ .* data.directions)
+        else
+            l_emp = norm(data.kappas)^2 - 3.0 * mean(data.kappas .* u_ .* data.directions)
+        end
+        # Regularization
+        l_reg = 0.0
+        if !isnothing(params.reg)
+            # for (order, power, λ) in params.reg    
+            for reg in params.reg    
+                # l_reg += reg.λ * regularization(θ; order=reg.order, power=reg.power)      
+                l_reg += regularization(θ, reg)      
+            end 
+        end
+        return l_emp + l_reg
+    end
+
+    function regularization(θ::ComponentVector, reg::AbstractRegularization; n_nodes=100)
+        
+        # Create (uniform) spacing time
+        # Δt = (params.tmax - params.tmin) / n_nodes
+        # times_reg = collect(params.tmin:Δt:params.tmax)
+
+        # LinRange does not propagate thought the backprop step!
+        # times_reg = collect(LinRange(params.tmin, params.tmax, n_nodes))
+        l_ = 0.0
+        if reg.order==0
+            l_ += quadrature(t -> norm(U([t], θ, st)[1])^reg.power, params.tmin, params.tmax, n_nodes)
+            # for t in times_reg
+            #     l_ += norm(U([t], θ, st)[1])^reg.power
+            # end
+        elseif reg.order==1
+            if reg.diff_mode=="AD"
+                # Compute gradient using automatic differentiaion in the NN
+                # @infiltrate
+                l_ += quadrature(t -> norm(Array(batched_jacobian(AutoForwardDiff(), τ -> U([τ], θ, st)[1],[t])))^reg.power, params.tmin, params.tmax, n_nodes)
+                # This currently doesn't run... too slow.
+                # for t in times_reg
+                #     # Try ReverseDiff
+                #     grad = Zygote.jacobian(first ∘ U, [t], θ, st)[1]
+                #     l_ += norm(grad)^reg.power
+                # end
+            elseif reg.diff_mode=="Finite Differences"
+                # Compute finite differences
+                L_estimated = map(t -> (first ∘ U)([t], θ, st), times_reg)
+                dLdt = diff(L_estimated) ./ diff(times_reg)
+                l_ += sum(norm.(dLdt).^reg.power)
+            else
+                throw("Method not implemented.")
+            end
+        else
+            throw("Method not implemented.")
+        end
+        return reg.λ * l_
+    end
+
+    losses = Float64[]
+    callback = function (p, l)
+        push!(losses, l)
+        if length(losses) % 200 == 0
+            println("Current loss after $(length(losses)) iterations: $(losses[end])")
+        end
+        return false
+    end
+
+    # adtype = Optimization.AutoZygote()
+    adtype = Optimization.AutoForwardDiff()
+    optf = Optimization.OptimizationFunction((x, θ) -> loss(x), adtype)
+    optprob = Optimization.OptimizationProblem(optf, ComponentVector{Float64}(θ))
+
+    res1 = Optimization.solve(optprob, ADAM(0.002), callback=callback, maxiters=params.niter_ADAM)
+    println("Training loss after $(length(losses)) iterations: $(losses[end])")
+
+    optprob2 = Optimization.OptimizationProblem(optf, res1.u)
+    res2 = Optimization.solve(optprob2, Optim.LBFGS(), callback=callback, maxiters=params.niter_LBFGS)
+    println("Final training loss after $(length(losses)) iterations: $(losses[end])")
+
+    # Optimized NN parameters
+    θ_trained = res2.u
+
+    # Final Fit 
+    fit_times = collect(params.tmin:0.1:params.tmax)
+    fit_directions = predict(θ_trained, T=fit_times)
+
+    return Results(θ_trained=θ_trained, U=U, st=st,
+                   fit_times=fit_times, fit_directions=fit_directions)
+end
diff --git a/test_complex_diff.jl b/test_complex_diff.jl
new file mode 100644
index 0000000..f044d4f
--- /dev/null
+++ b/test_complex_diff.jl
@@ -0,0 +1,25 @@
+using Lux, Zygote
+using Random
+rng = Random.default_rng()
+Random.seed!(rng, 666)
+
+rbf(x) = exp.(-(x .^ 2))
+
+U = Lux.Chain(
+    Lux.Dense(1, 10, rbf),
+    Lux.Dense(10, 3, rbf)
+)
+
+θ, st = Lux.setup(rng, U)
+
+@show U([1.0], θ, st)[begin]
+
+function complex_step_differentiation(f::Function, x::Float64, ϵ::Float64)
+    return imag(f(x + ϵ * im)) / ϵ
+end
+
+@show complex_step_differentiation(t -> U([t], θ, st)[begin], 1.0, 1e-5)
+
+loss(t) = complex_step_differentiation(τ -> U([τ], θ, st)[begin], t, 1e-5)
+
+Zygote.gradient(t -> loss(t), 1.0)
\ No newline at end of file
diff --git a/test_serialization.jl b/test_serialization.jl
new file mode 100644
index 0000000..7a8c05a
--- /dev/null
+++ b/test_serialization.jl
@@ -0,0 +1,25 @@
+using JLD2
+
+abstract type AbstractData end
+
+@kwdef struct Data{F <: AbstractFloat} <: AbstractData
+    times::Vector{F}
+    directions::Matrix{F}
+    kappas::Union{Vector{F}, Nothing}
+    L::Union{Function, Nothing}
+end
+
+data = Data(times=ones(4), directions=ones(4,2), kappas=nothing, L= x-> x^3)
+
+@kwdef struct DataSerialization{F <: AbstractFloat}
+    times::Vector{F}
+    directions::Matrix{F}
+    kappas::Union{Vector{F}, Nothing}
+    L::Union{Function, Nothing}
+end
+
+JLD2.writeas(::Type{Data}) = DataSerialization
+Base.convert(::Type{DataSerialization}, a::Data) = DataSerialization([a.times, a.directions, a.kappas, a.L])
+Base.convert(::Type{Data}, a::DataSerialization) = Data(times = a.times, directions = a.directions, kappas=a.kappas, L=a.L)
+
+@save "test_serialization.jld2" data
\ No newline at end of file
diff --git a/test_serialization.jld2 b/test_serialization.jld2
new file mode 100644
index 0000000..f906b8c
Binary files /dev/null and b/test_serialization.jld2 differ