#Packages using LinearAlgebra using Plots using FFTW using Printf using BenchmarkTools using Rotations using StaticArrays #%Complex definition i=sqrt(Complex(-1)); #Definition of individual spin operators nspins=1 include("basis.jl") iu,ix,iy,iz,ip,im,ia,ib=basis(nspins) #General Setup timestep = 0.000001 w = StaticVector{3,Float64} w = [1.0, 1000, 0.0] # z can't be zero! eps = 0.000000001 th = w * timestep ham = 2 * pi * th[1] * ix[:,:,1] + 2 * pi * th[2] * iy[:,:,1] + 2 * pi * th[3] * iz[:,:,1] uham = exp(- i * ham) rho = iy[:,:,1] rho = uham * rho * uham' function derivative_x_expm(w::Vector{Float64},timestep::Float64) local ham = Array{ComplexF64,2}(undef,2,2) local ctrl= Array{ComplexF64,2}(undef,2,2) local res = Array{ComplexF64,2}(undef,4,4) local i = sqrt(ComplexF64(-1)); local ham = 2 * pi * w[1] * ix[:,:,1] + 2 * pi * w[2] * iy[:,:,1] + 2 * pi * w[3] * iz[:,:,1] local ctrl = - 2 * pi * ix[:,:,1] ./(-timestep*2*pi) local mnull = zeros(2,2) res = exp(-i * timestep * [ham ctrl; mnull ham]) return res end function derivative_x_xact(w::Vector{Float64},timestep::Float64) dQdx = Vector{Float64}(undef,4) x = Float64 y = Float64 z = Float64 r = Float64 si = Float64 co = Float64 local x,y,z = w[1:3] .* (2 * pi * timestep ) local r = sqrt(x^2 + y^2 + z^2) local si,co = sincos(r) local dRdx = zeros(3, 3) dRdx[1, 1] = (x^3 / r^3 - x / r) * si - (x / r^2 - x^3 / r^4) * 2 * (co - 1) dRdx[2, 2] = (x * y^2 / r^3 - x / r) * si + (x * y^2 / r^4) * 2 * (co - 1) dRdx[3, 3] = (x * z^2 / r^3 - x / r) * si + (x * z^2 / r^4) * 2 * (co - 1) dRdx[1, 2] = (x^2 * y + x * z) / r^3 * si + (2 * x^2 * y / r^4 - y / r^2) * (co - 1) - (x * z / r^2) * co dRdx[2, 1] = (x^2 * y - x * z) / r^3 * si + (2 * x^2 * y / r^4 - y / r^2) * (co - 1) + (x * z / r^2) * co dRdx[1, 3] = (x^2 * y - x * y) / r^3 * si + (2 * x^2 * z / r^4 - y / r^2) * (co - 1) + (x * y / r^2) * co dRdx[3, 1] = (x^2 * y + x * y) / r^3 * si + (2 * x^2 * z / r^4 - y / r^2) * (co - 1) - (x * y / r^2) * co dRdx[2, 3] = (x * y * z - z^2 - y^2) / r^3 * si + 2 * x * y * z / r^4 * (co - 1) - x^2 / r^2 * co dRdx[3, 2] = (x * y * z + z^2 + y^2) / r^3 * si + 2 * x * y * z / r^4 * (co - 1) + x^2 / r^2 * co return dRdx end function derivative_alpha_xact(w::Vector{Float64},timestep::Float64) dRdalpha = Array{Float64,2}(undef,3,3) x = Float64 y = Float64 z = Float64 r = Float64 si = Float64 co = Float64 sinalpha = Float64 cosalpha = Float64 local txy, z = w[2:3] .* (2 * pi * timestep ) local r = sqrt(txy^2 + z^2) local si,co = sincos(r) local sinalpha, cosalpha = sincos(w[1]) local nx = txy * cosalpha/r local ny = txy * sinalpha/r local nz = z/r local dRdalpha = zeros(3, 3) dRdalpha[1,1] = -2 * nx * ny * (1 - co) dRdalpha[2,2] = 2 * nx * ny * (1 - co) dRdalpha[3,3] = 0 dRdalpha[1,2] = (nx^2 - ny^2) * (1 - co) dRdalpha[2,1] = (nx^2 - ny^2) * (1 - co) dRdalpha[1,3] = -ny * nz * (1 - co) + nx * si dRdalpha[3,1] = -ny * nz * (1 - co) - nx * si dRdalpha[2,3] = nx * nz * (1 - co) + ny * si dRdalpha[3,2] = nx * nz * (1 - co) - ny * si return dRdalpha end function derivative_thetaxy_xact(w::Vector{Float64},timestep::Float64) dRdtxy = Array{Float64,2}(undef,3,3) x = Float64 y = Float64 z = Float64 r = Float64 si = Float64 co = Float64 sinalpha = Float64 cosalpha = Float64 nx = Float64 ny = Float64 nz = Float64 local txy, z = w[2:3] .* (2 * pi * timestep ) local r = sqrt(txy^2 + z^2) local si,co = sincos(r) local sinalpha, cosalpha = sincos(w[1]) local nx = txy * cosalpha/r local ny = txy * sinalpha/r local nz = z/r local dRdtxy = zeros(3, 3) dRdtxy[1,1] = -txy * (ny^2 + nz^2) / r * si + 2 * nx^2 * nz^2 / txy * (1 - co) dRdtxy[2,2] = -txy * (nx^2 + nz^2) / r * si + 2 * ny^2 * nz^2 / txy * (1 - co) dRdtxy[3,3] = -txy * nz^2 / r^3 * (txy^2 * si + 2 * (1 + co)) - txy^5 / r^5 * si dRdtxy[1,2] = txy * nz / r^2 * (si - r * co) + nx * ny / (txy * r^2) * (txy^2 * r * si + (1 - co) * (r^2 - txy^2)) dRdtxy[2,1] = txy * nz / r^2 * (r * co - si) + nx * ny / (txy * r^2) * (txy^2 * r * si + (1 - co) * (r^2 - txy^2)) dRdtxy[1,3] = ny / (txy * r^2) * ((r^2 - txy^2) * si + txy^2 * r * co) + nx / (txy * r^3) * ((r^2 * z - 2 * z * txy^2) * (1 - co) + txy^2 * z * si) dRdtxy[3,1] = ny / (txy * r^2) * ((txy^2 - r^2) * si - txy^2 * r * co) + nx / (txy * r^3) * ((r^2 * z - 2 * z * txy^2) * (1 - co) + txy^2 * z * si) dRdtxy[2,3] = nx / (txy * r^2) * ((txy^2 - r^2) * si - txy^2 * r * co) + ny / (txy * r^3) * ((r^2 * z - 2 * z * txy^2) * (1 - co) + txy^2 * z * si) dRdtxy[3,2] = nx / (txy * r^2) * ((r^2 - txy^2) * si + txy^2 * r * co) + ny / (txy * r^3) * ((r^2 * z - 2 * z * txy^2) * (1 - co) + txy^2 * z * si) return dRdtxy end function derivative_z_xact(w::Vector{Float64},timestep::Float64) dRdz = Array{Float64,2}(undef,3,3) r = Float64 si = Float64 co = Float64 sinalpha = Float64 cosalpha = Float64 nx = Float64 ny = Float64 nz = Float64 local txy, z = w[2:3] .* (2 * pi * timestep ) local r = sqrt(txy^2 + z^2) local si,co = sincos(r) local sinalpha, cosalpha = sincos(w[1]) local nx = txy * cosalpha/r local ny = txy * sinalpha/r local nz = z/r local dRdz = zeros(3, 3) dRdz[1,1] = (nx^2 * nz - (nz) ) * si + 2 * nx^2 * nz * (co - 1) / r dRdz[1,2] = (nx * ny * nz - txy^2 / r^3) * si + 2 * nx * ny * nz / r * (co - 1) - nz^2 * co dRdz[1,3] = (nx * nz^2 - ny * nz / r) * si + (2 * nx * nz^2 - nx) * (co - 1) / r + ny * nz * co dRdz[2,1] = (nx * ny * nz + txy^2 / r^3) * si + 2 * nx * ny * nz / r * (co - 1) + nz^2 * co dRdz[2,2] = (ny^2 * nz - nz) * si - 2 * ny^2 * nz * (1 - co) / r dRdz[2,3] = (ny * nz^2 + nx * nz / r) * si + (2 * ny * nz^2 - ny) * (co - 1) / r - nx * nz * co dRdz[3,1] = (nx * nz^2 + ny * nz / r) * si + (2 * nx * nz^2 - nx) * (co - 1) / r - ny * nz * co dRdz[3,2] = (ny * nz^2 - nx * nz / r) * si + (2 * ny * nz^2 - ny) * (co - 1) / r + nx * nz * co dRdz[3,3] = (nz^3 - nz) * si + (nz - nz^3) * 2 * (1 - co) / r return dRdz end function derivative_alpha_diff(w::Vector{Float64},timestep::Float64,eps::Float64) R = Array{Float64,2}(undef,3,3) Re = Array{Float64,2}(undef,3,3) dR = Array{Float64,2}(undef,3,3) r = Float64 alpha = Float64 txy = Float64 si = Float64 co = Float64 sinalpha = Float64 cosalpha = Float64 nx = Float64 ny = Float64 nz = Float64 local txy, z = w[2:3] .* (2 * pi * timestep ) local alpha = w[1] local r = sqrt(txy^2 + z^2) local si,co = sincos(r) local sinalpha, cosalpha = sincos(alpha) local nx = txy * cosalpha/r local ny = txy * sinalpha/r local nz = z/r local R = zeros(3,3) R[1,1] = co + nx^2 * (1 - co) R[1,2] = - nz * si + nx * ny * (1 - co) R[1,3] = ny * si + nx * nz * (1 - co) R[2,1] = nz * si + nx * ny * (1 - co) R[2,2] = co + ny^2 * (1 - co) R[2,3] = - nx * si + ny * nz * (1 - co) R[3,1] = -ny * si + nx * nz * (1 - co) R[3,2] = nx * si + ny * nz * (1 - co) R[3,3] = co + nz^2 * (1 - co) alpha = alpha + eps local sinalpha, cosalpha = sincos(alpha) local nx = txy * cosalpha/r local ny = txy * sinalpha/r # local si,co = sincos(r) local Re = zeros(3,3) Re[1,1] = co + nx^2 * (1 - co) Re[1,2] = - nz * si + nx * ny * (1 - co) Re[1,3] = ny * si + nx * nz * (1 - co) Re[2,1] = nz * si + nx * ny * (1 - co) Re[2,2] = co + ny^2 * (1 - co) Re[2,3] = - nx * si + ny * nz * (1 - co) Re[3,1] = -ny * si + nx * nz * (1 - co) Re[3,2] = nx * si + ny * nz * (1 - co) Re[3,3] = co + nz^2 * (1 - co) local dR = (Re-R)/eps return dR end function derivative_thetaxy_diff(w::Vector{Float64},timestep::Float64,eps::Float64) R = Array{Float64,2}(undef,3,3) Re = Array{Float64,2}(undef,3,3) dR = Array{Float64,2}(undef,3,3) r = Float64 alpha = Float64 txy = Float64 si = Float64 co = Float64 sinalpha = Float64 cosalpha = Float64 nx = Float64 ny = Float64 nz = Float64 local txy, z = w[2:3] .* (2 * pi * timestep ) local alpha = w[1] local r = sqrt(txy^2 + z^2) local si,co = sincos(r) local sinalpha, cosalpha = sincos(alpha) local nx = txy * cosalpha/r local ny = txy * sinalpha/r local nz = z/r local R = zeros(3,3) R[1,1] = co + nx^2 * (1 - co) R[1,2] = - nz * si + nx * ny * (1 - co) R[1,3] = ny * si + nx * nz * (1 - co) R[2,1] = nz * si + nx * ny * (1 - co) R[2,2] = co + ny^2 * (1 - co) R[2,3] = - nx * si + ny * nz * (1 - co) R[3,1] = -ny * si + nx * nz * (1 - co) R[3,2] = nx * si + ny * nz * (1 - co) R[3,3] = co + nz^2 * (1 - co) txy = txy + eps local r = sqrt(txy^2 + z^2) local si,co = sincos(r) local nx = txy * cosalpha/r local ny = txy * sinalpha/r local nz = z/r local Re = zeros(3,3) Re[1,1] = co + nx^2 * (1 - co) Re[1,2] = - nz * si + nx * ny * (1 - co) Re[1,3] = ny * si + nx * nz * (1 - co) Re[2,1] = nz * si + nx * ny * (1 - co) Re[2,2] = co + ny^2 * (1 - co) Re[2,3] = - nx * si + ny * nz * (1 - co) Re[3,1] = -ny * si + nx * nz * (1 - co) Re[3,2] = nx * si + ny * nz * (1 - co) Re[3,3] = co + nz^2 * (1 - co) local dR = (Re-R)/eps return dR end function derivative_z_diff(w::Vector{Float64},timestep::Float64,eps::Float64) R = Array{Float64,2}(undef,3,3) Re = Array{Float64,2}(undef,3,3) dR = Array{Float64,2}(undef,3,3) r = Float64 alpha = Float64 txy = Float64 si = Float64 co = Float64 sinalpha = Float64 cosalpha = Float64 nx = Float64 ny = Float64 nz = Float64 local txy, z = w[2:3] .* (2 * pi * timestep ) local alpha = w[1] local r = sqrt(txy^2 + z^2) local si,co = sincos(r) local sinalpha, cosalpha = sincos(alpha) local nx = txy * cosalpha/r local ny = txy * sinalpha/r local nz = z/r local R = zeros(3,3) R[1,1] = co + nx^2 * (1 - co) R[1,2] = - nz * si + nx * ny * (1 - co) R[1,3] = ny * si + nx * nz * (1 - co) R[2,1] = nz * si + nx * ny * (1 - co) R[2,2] = co + ny^2 * (1 - co) R[2,3] = - nx * si + ny * nz * (1 - co) R[3,1] = -ny * si + nx * nz * (1 - co) R[3,2] = nx * si + ny * nz * (1 - co) R[3,3] = co + nz^2 * (1 - co) z = z + eps local r = sqrt(txy^2 + z^2) local nx = txy * cosalpha/r local ny = txy * sinalpha/r local nz = z/r local Re = zeros(3,3) Re[1,1] = co + nx^2 * (1 - co) Re[1,2] = - nz * si + nx * ny * (1 - co) Re[1,3] = ny * si + nx * nz * (1 - co) Re[2,1] = nz * si + nx * ny * (1 - co) Re[2,2] = co + ny^2 * (1 - co) Re[2,3] = - nx * si + ny * nz * (1 - co) Re[3,1] = -ny * si + nx * nz * (1 - co) Re[3,2] = nx * si + ny * nz * (1 - co) Re[3,3] = co + nz^2 * (1 - co) local dR = (Re-R)/eps return dR end function printcmat(mat) nrow,ncol = size(mat) for irow = 1:nrow @printf(" ") for icol = 1:ncol @printf("(%f %f*i) ",real(mat[irow,icol]), imag(mat[irow,icol])) end @printf("\n") end end testxact1 = derivative_alpha_xact(w,timestep) testdiff1 = derivative_alpha_diff(w,timestep,eps) testxact2 = derivative_thetaxy_xact(w,timestep) testdiff2 = derivative_thetaxy_diff(w,timestep,eps) testxact3 = derivative_z_xact(w,timestep) testdiff3 = derivative_z_diff(w,timestep,eps) println("dRdalpha xact: ") printcmat(testxact1[1:3,1:3]) println("dRdalpha diff: ") printcmat(testdiff1[1:3,1:3]) println("dRdxy xact: ") printcmat(testxact2[1:3,1:3]) println("dRdxy diff: ") printcmat(testdiff2[1:3,1:3]) println("dRdz xact: ") printcmat(testxact3[1:3,1:3]) println("dRdz diff: ") printcmat(testdiff3[1:3,1:3]) ntest = 1000 wtest = rand(Float64,ntest,3) println("Benchmark for xact solution alpha,txy,z: ") @btime begin for itest = 1:ntest testxact1 = derivative_alpha_xact(wtest[itest,:],timestep) end end @btime begin for itest = 1:ntest testxact2 = derivative_thetaxy_xact(wtest[itest,:],timestep) end end @btime begin for itest = 1:ntest testxact2 = derivative_z_xact(wtest[itest,:],timestep) end end println("Benchmark for finite difference solution alpha,txy,z: ") @btime begin for itest = 1:ntest testdiff1 = derivative_alpha_diff(wtest[itest,:],timestep,eps) end end @btime begin for itest = 1:ntest testdiff1 = derivative_thetaxy_diff(wtest[itest,:],timestep,eps) end end @btime begin for itest = 1:ntest testdiff1 = derivative_z_diff(wtest[itest,:],timestep,eps) end end