Common Propagator types

common_type.jl

using QuBase
using Docile

import Base: start, next, done

abstract QuPropagatorMethod

immutable QuPropagator{QPM<:QuPropagatorMethod}
    hamiltonian
    init_state::QuVector
    tlist
    method::QPM
    QuPropagator(hamiltonian, init_state, tlist, method) = new(hamiltonian, init_state, tlist, method)
end

QuPropagator{QPM<:QuPropagatorMethod}(hamiltonian, init_state, tlist, method::QPM) = QuPropagator{QPM}(hamiltonian, init_state, tlist, method)

immutable QuPropagatorState
    psi::QuVector
    t
    t_state
end

function Base.start(prob::QuPropagator)
    t_state = start(prob.tlist)
    t,t_state = next(prob.tlist,t_state)
    return QuPropagatorState(prob.init_state,t,t_state)
end

Base.done(prob::QuPropagator, qustate::QuPropagatorState) = done(prob.tlist, qustate.t_state)

function Base.next{QPM<:QuPropagatorMethod}(prob::QuPropagator{QPM}, qustate::QuPropagatorState)
    op = prob.hamiltonian
    current_qustate = qustate.psi
    current_t = qustate.t
    if done(prob, qustate) == false
        t,t_state = next(prob.tlist, qustate.t_state)
        next_qustate = propagate(prob.method, op, t, current_t, current_qustate)
        return (next_qustate, t), QuPropagatorState(next_qustate, t, t_state)
    else
        return QuPropagatorState(current_qustate, current_t, qustate.t_state)
    end
end

immutable QuEuler <: QuPropagatorMethod
end

function propagate(prob::QuEuler, op, t, current_t, current_qustate)
    dt = t - current_t
    return (eye(op)-im*op*dt)*current_qustate
end

immutable QuCN <: QuPropagatorMethod
end

function propagate(prob::QuCN, op, t, current_t, current_qustate)
    dt = t - current_t
    uni = eye(op)-im*op*dt/2
    return \(uni', uni*current_qustate)
end

immutable QuKrylov <: QuPropagatorMethod
    options::Dict
end

function propagate(prob::QuKrylov, op, t, current_t, current_qustate)
    basis_size = prob.options[:NC]
    dt = t-current_t
    N = min(basis_size, length(current_qustate))
    v = Array(typeof(current_qustate),0)
    push!(v,zeros(current_qustate))
    push!(v,current_qustate)
    w = Array(typeof(current_qustate),0)
    alpha = Array(Complex{Float64},0)
    beta = Array(Complex{Float64},0)
    push!(beta,0.)
    for i=2:N
        push!(w, op*v[i])
        push!(alpha, w[i-1]'*v[i])
        w[i-1] = w[i-1]-alpha[i-1]*v[i]-beta[i-1]*v[i-1]
        push!(beta, norm(w[i-1]))
        push!(v, w[i-1]/beta[i])
    end
    push!(w, op*v[end])
    push!(alpha, w[end]'*v[end])
    deleteat!(v,1)
    H_k = full(Tridiagonal(beta[2:end], alpha, beta[2:end]))
    U_k = expm(-im*dt*H_k)
    next_state = zeros(current_qustate)
    for j=1:N
        next_state = next_state + v[j]*U_k[j,1]
    end
    return next_state
end