How to plot a time-evolution of 3D Gaussian in Julia?

Viewed 228

I want to plot a time-evolution of 3D Gaussian with Makie.jl. Here is a surface-version code of sin(r)/r.
So I wrote a code in reference to it.

using Makie  
using FileIO  
using LinearAlgebra  
using AbstractPlotting

scene = Scene(backgroundcolor = :black);  
f(x,y,z) = exp(-((x)^2 + (y)^2 + (z)^2))  
r = LinRange(-5, 5, 50)  
vol_func(t) = [Float64(f(x - cos(t),y - sin(t),z - t)) for x = r, y = r,z = r]  
vol = volume!(scene,r,r,r,vol_func(20),algorithm = :mip)[end]  
scene[Axis].names.textcolor = :gray  
N = 20  
scene 
record(scene, "voloutput.mp4", range(0, stop = 5, length = N)) do t  
    vol[3] = vol_func(t)  
end  

But this code does not work.

MethodError: Cannot `convert` an object of type Array{Float64,3} to an object of type LinRange{Float64}

How should I fix the code?

P.S.
The snapshot at initial time is like this.(reference)

using Makie
using FileIO
using LinearAlgebra
using AbstractPlotting

r = LinRange(-20, 20, 500);  # our value range

ρ(x, y, z) = exp(-((x-1)^2 + (y)^2 + (z)^2)) # function (charge density)

# create a Scene with the attribute `backgroundcolor = :black`,
# can be any compatible color.  Useful for better contrast and not killing your eyes with a white background.
scene = Scene(backgroundcolor = :black)

volume!(
    scene,
    r, r, r,          # coordinates to plot on
    ρ,                # charge density (functions as colorant)
    algorithm = :mip  # maximum-intensity-projection
)

scene[Axis].names.textcolor = :gray # let axis labels be seen on dark 
background

save("sp.png",scene)

I want to see the yellow region moving as spiral. (2020/08/28)

I just realized not vol[3] but vol[4]. Then, it worked. But I have a next question. (2020/08/31)
I tried to do the same thing for the matrix-form time-dependent Schrodinger equation with its initial condition being Gaussian.

using LinearAlgebra
using OrdinaryDiffEq
using DifferentialEquations  
  
#Define the underlying equation
function time_evolution(ψdot,ψ,p,t)
  ψdot.=-im.*H(Lx,Ly,Lz)*ψ
end

Lx = Ly = Lz = 10
ψ0 = [] #  Initial conditions


for iz = 1:Lz
    for ix = 1:Lx
        for iy = 1:Ly                  
           gauss = exp(-((ix)^2 + (iy)^2 + (iz)^2))
           push!(ψ0,gauss)                           
        end
    end
end

tspan = (0.,1.0) #  Simulation time span


#Pass to Solvers
prob = ODEProblem(time_evolution,ψ0,tspan)
sol = solve(prob)

Here,H(Lx,Ly,Lz) is a N×N matrix parameterized by systemsize Lx,Ly,Lz and N = Lx×Ly×Lz. The sample code of H(Lx,Ly,Lz) is here. Then,

using Makie
using FileIO 
using LinearAlgebra
using AbstractPlotting
using ColorSchemes

x = 1: Lx  # our value range
y = 1: Ly
z = 1: Lz

ρ(ix,iy,iz,nt) = abs2.((sol[nt][(iz-1)*Lx*Ly + (ix-1)*Ly + (iy-1)])./norm(sol[nt][(iz-1)*Lx*Ly + (ix-1)*Ly + (iy-1)]))
ψ(nt) = Float64[ρ(ix,iy,iz,nt) for ix in x, iy in y,iz in z]
scene = Scene(backgroundcolor = :white)

c = ψ(length(sol.t))

vol = volume!(
scene,
x, y, z,          # coordinates to plot on
c,                # charge density (functions as colorant)
algorithm = :mip, # maximum-intensity-projection
colorrange = (0,0.01),
transparency = true,
)[end]

update_cam!(scene, Vec3f0(1,0.5,0.1), Vec3f0(0))
scene[Axis].names.textcolor = :gray # let axis labels be seen on darkbackground

record(scene, "output.mp4", range(0, stop = length(sol.t)-1, length = 1)) do nt
    vol[4] = ψ(nt)
end

But this code has an error.

ArgumentError: range(0.0, stop=5.0, length=1): endpoints differ

Where is the mistake?

I found the mistake.(2020/09/02)

sol[nt]→sol(nt)
range(0, stop = length(sol.t)-1, length = 1)→range(0, stop = 1.0, length = 20)

Then, the code passed and a mp4 animation was obtained.
But the plot can't be seen in the mp4 file. Why...

0 Answers
Related