--- a/lattice2d.jl	Sun Jan 09 11:27:04 2022 +0100
+++ b/lattice2d.jl	Wed Jan 19 22:51:26 2022 +0100
@@ -4,7 +4,7 @@
 
 using SparseArrays
 
-theme(:juno)
+theme(:ggplot2)
 
 # Masses schema:
 #
@@ -82,25 +82,25 @@
 
         # Next-neighbor out-dimensional coupling: i.e., top neighbor causes x displacement.
         # The coupling constant is just half here. In reality, there would be a nonlinear dependency.
-        A[2N+x, 2t-1] += k/4
-        A[2N+x, 2b-1] += k/4
-        A[2N+x, 2i-1] -= k/2
+        # A[2N+x, 2t-1] += k/4
+        # A[2N+x, 2b-1] += k/4
+        # A[2N+x, 2i-1] -= k/2
 
-        A[2N+y, 2r] += k/4
-        A[2N+y, 2l] += k/4
-        A[2N+y, 2i] -= k/2
+        # A[2N+y, 2r] += k/4
+        # A[2N+y, 2l] += k/4
+        # A[2N+y, 2i] -= k/2
 
         # Add diagonal neighbors
-        #
-        # for e in (r1, r2, l1, l2)
-        #     # Diagonal dependence on √2 neighbors.
-        #     A[2N+x, 2*e-1] += k/4
-        #     A[2N+x, 2*e]   += k/8 # Interaction between diag. neighbor y and this element's x is weaker.
-        #     A[2N+y, 2*e-1] += k/8
-        #     A[2N+y, 2*e]   += k/4
-        # end
-        # A[2N+x, x] -= 3/2 * k
-        # A[2N+y, y] -= 3/2 * k
+        
+        for e in (r1, r2, l1, l2)
+            # Diagonal dependence on √2 neighbors.
+            A[2N+x, 2*e-1] += k/4
+            A[2N+x, 2*e]   += k/8 # Interaction between diag. neighbor y and this element's x is weaker.
+            A[2N+y, 2*e-1] += k/8
+            A[2N+y, 2*e]   += k/4
+        end
+        A[2N+x, x] -= 3/2 * k
+        A[2N+y, y] -= 3/2 * k
 
         # Add damping terms.
         A[2N+1:end, 2N+1:end] -= LinearAlgebra.I(2N) * damping
@@ -112,9 +112,8 @@
     u0 = zeros(4n^2)
     N = n^2
     for (x, y, dx, dy) in displacements
-        i = n*(x-1) + y
+        i = n*(y-1) + x
         ix, iy = 2i-1, 2i
-        #@show (i, ix, iy)
         u0[ix] += dx
         u0[iy] += dy
     end
@@ -149,7 +148,7 @@
     [(i, 2, .5, .0) for i in 1:n])
 initial_xy_transversal_displacements = vcat(
     [(1, i, 0.5, 0.0) for i in 1:n],
-    [(i, 2, 0.0, 0.5) for i in 1:n])
+    [(i, 1, 0.0, 0.5) for i in 1:n])
 initial_x_longitudinal_displacements = [(i, 1, .2, 0.) for i in 1:n]
 initial_x_transversal_displacements = [(i, 1, 0., 0.2) for i in 1:n]
 #  One element is displaced diagonally.
@@ -160,10 +159,14 @@
 u0 = initial_state(n, initial_xy_transversal_displacements)
 
 # Essential part: the coupling matrix!
-coupling = coupling_matrix(n, 1.; damping=0.0000)
+coupling = coupling_matrix(n, 1.; damping=0.0000, full=true)
+
+# Alternatively: Use an eigenvector as initial state, exciting an eigenmode.
+ED = LinearAlgebra.eigen(coupling[2n^2+1:end, 1:2n^2]);
+u1 = vcat(ED.vectors[:,150] ./ 2, zeros(2n^2));
 
 # ODE formalities
-problem = ODEProblem(f2d, u0, (0, 20), coupling);
+problem = ODEProblem(f2d, u0, (0, 50), coupling);
 @time solution = solve(problem);
 
 @gif for (j, u) in enumerate(solution.u)

--- a/linear_chain.jl	Sun Jan 09 11:27:04 2022 +0100
+++ b/linear_chain.jl	Wed Jan 19 22:51:26 2022 +0100
@@ -2,7 +2,7 @@
 import LinearAlgebra
 using DifferentialEquations
 using Plots
-theme(:juno)
+theme(:ggplot2)
 
 # The coupling matrix has 4nx4n elements for a linear chain of two degrees of freedom per oscillator.
 # We use Born-von Karman boundary conditions, i.e. the chain is a ring.
@@ -97,30 +97,48 @@
     gif(anim, "anim_$(n).gif")
 end
 
+function plot_website_chain_animation(n::Int, solution)
+    xs = .4 .* collect(1:n)
+    ys = zeros(n)
+    T = length(solution.u)
+
+    trajectories_x, trajectories_y = zeros(n, T), zeros(n, T)
+
+    anim = @animate for (i,u) in enumerate(solution.u)
+        xpos = u[1:2:2n]
+        ypos = u[2:2:2n]
+        trajectories_x[:, i] = xs .+ xpos
+        trajectories_y[:, i] = ypos .+ ones(n)
+        plot(vcat([0], xs .+ xpos, [(n+1)*.4]), ones(n+2) .+ vcat([0], ypos, [0]), xlim=(0, .4*(n+1)), ylim=(.8, 1.2), size=(800, 200), legend=false, markersize=8, markershape=:circle)
+        plot!(trajectories_x[:, 1:i]', trajectories_y[:, 1:i]')
+    end
+    gif(anim, "anim_$(n).gif", fps=12)
+end
+
 function initial_wave_displacement(N::Int, q::Float64, ampl::Float64=.25, a::Float64=1.; trans=false)::Vector{Tuple}
     N = div(N, 2)
     ampls = [ampl * sin(q * a * i) for i in 1:N]
     [(i, trans ? 0. : ampls[i], trans ? ampls[i] : 0.) for i in 1:N]
 end
 
-N = 15
-k1, k2 = .4, 1
-mass = 1
-initial_displacement = [(1, .0, .3), (N, -.0, .0)]
+N = 20
+k1, k2 = 1., 2.
+mass = 1.
+initial_displacement = [(1, 0.15, 0.15), (N, 0., 0.)]
 initial_displacement_2 = [(1, .5, 0.)]
 initial_wave_4 = initial_wave_displacement(N, 2pi/(N), .1; trans=true)
 timespan = (0, 90)
-coupling_H = n_coupled_masses(N, k1 * ones(N+1), mass * ones(N); fixed_ends=false)
+coupling_H = n_coupled_masses(N, k1 * ones(N+1), ones(N) * mass; fixed_ends=true)
 
 E = eigenmodes(N, coupling_H)
 
-# Two different initial configurations as choice.
+# Three different initial configurations as choice.
 u0 = initial_state(N, initial_displacement)
 u1 = initial_state_from_eigenmode(E.vectors[:,1] ./ 3)
 u2 = initial_state(N, initial_wave_4)
 
-problem = ODEProblem(f, u0, timespan, coupling_H)
-@time solution = solve(problem)
+problem = ODEProblem(f, u0, timespan, coupling_H);
+@time solution = solve(problem);
 
 # 3D trajectory plot
 # plot(solution, vars=[(0, 2x-1, 2x) for x in 1:N],