Till 27-05-2015

270515.jl

using QuBase
using Docile

import Base: start, next, done

abstract QuantumPropagator{B<:AbstractBasis,T,A}

type StepPropagator{B<:AbstractBasis,T,A} <: QuantumPropagator{B,T,A}
    hamiltonian::QuMatrix{B,T,A}
    init_state::QuVector
    tlist::Array
    StepPropagator(hamiltonian,init_state,tlist) = new(hamiltonian,init_state,tlist)
end

type StepPropagatorState
    psi::QuVector
    t_index::Int
end

immutable TimeStepState
    tlist::Array
    t_index::Int
end

Base.start(state::TimeStepState) = state.tlist[1]

Base.next(state::TimeStepState) = state.tlist[state.t_index+1], state.tlist[state.t_index]

Base.done(state::TimeStepState) = state.tlist[state.t_index] == state.tlist[end] ? true : false

@doc """
Euler Problem
`hamiltonian` : Hamiltonian of the system
`init_state`  : initial state of the syste
`tlist`       : Time step array
""" ->
immutable QuEuler{B<:AbstractBasis,T,A} <: QuantumPropagator{B,T,A}
    hamiltonian::QuMatrix{B,T,A}
    init_state::QuVector
    tlist::Array
    QuEuler(hamiltonian, init_state, tlist) = new(hamiltonian, init_state, tlist)
end

QuEuler{B<:AbstractBasis,T,A}(hamiltonian::QuMatrix{B,T,A}, init_state::QuVector, tlist::Array)=QuEuler{B,T,A}(hamiltonian, init_state, tlist)
# Euler(hamiltonian, init_state, tlist) = solver_euler(hamiltonian, init_state, tlist)

@doc """
Euler State
`psi`     : State of the system
`t_index` : Array index of time in `tlist` of the current state of the system
""" ->
immutable QuEulerState
    psi::QuVector
    # t_index::Int
    t_state::TimeStepState
end

@doc """
Input Parameters : Euler Problem

Returns the starting, i.e., the initial state of the system
""" ->
function Base.start(prob::QuEuler)
    return QuEulerState(prob.init_state, TimeStepState(prob.tlist,1))
end

@doc """
Input Parameters : Euler Problem and Euler State

Returns the next state using Euler's method.
""" ->
function Base.next(prob::QuEuler, state::QuEulerState)
    op = prob.hamiltonian
    t_list = prob.tlist
    current_state = state.psi
    #current_t = state.t_index
    current_tstep = state.t_state
    #if t_list[current_t] != t_list[end]
    if current_tstep.tlist[current_tstep.t_index] != current_tstep.tlist[end]
        # next_t = t_list[current_t+1]
        #next_state = (eye(op)-im*op*(next_t-t_list[current_t]))*current_state
        #return (next_state, next_t), QuEulerState(next_state, findfirst(t_list, next_t))
        # next_t = next(current_step)[1]
        next_state = (eye(op)-im*op*(next(current_tstep)[1]-next(current_tstep)[2]))*current_state
        return (next_state, next(current_tstep)[1]), QuEulerState(next_state, TimeStepState(t_list,current_tstep.t_index+1))
    else
        # return QuEulerState(current_state, findfirst(t_list,current_t))
        return QuEulerState(current_state, TimeStepState(t_list,current_step.tindex))
    end
end

@doc """
Input Parameters : Euler Problem and Euler State

Returns true if the current state is final state, else false
""" ->
function Base.done(prob::QuEuler, state::QuEulerState)
    #if prob.tlist[state.t_index] == prob.tlist[end]
    # if done(state.t_state)
    #   return true
    # else
    #   return false
    #end
    return done(state.t_state) == true ? true : false
end


@doc """
Internal function to initialize the system
""" ->
function initialize(hamiltonian::QuMatrix, init_state::QuVector, tlist::Array)
    final_state = Array(typeof(init_state),0)
    sizehint(final_state, length(tlist))
    push!(final_state, init_state)
    return final_state
end

@doc """
Input Parameters : Hamiltonian of the system, Initial State of the system, Time steps

Returns the collection of output states at every time step using the Forward Euler Method.
""" ->
function solver_euler(hamiltonian::QuMatrix, init_state::QuVector, tlist::Array)
    final_state = initialize(hamiltonian, init_state, tlist)
    #for i=2:length(tlist)
    #    uni = QuArray(eye(size(hamiltonian,1)))-im*hamiltonian*(tlist[i]-tlist[i-1])
    #    push!(final_state, uni*final_state[i-1])
    #end
    id = eye(hamiltonian)
    next_step = start(tlist)
    while !done(tlist, next_step)
        step, next_step = next(tlist, next_step)
        if next_step == length(tlist)+1
          break
        else
          uni = id-im*hamiltonian*(tlist[next_step]-step)
          push!(final_state, uni*final_state[next_step-1])
        end
    end
    return final_state
end

@doc """
Crank Nicolson Problem
`hamiltonian` : Hamiltonian of the system
`init_state`  : initial state of the syste
`tlist`       : Time step array
""" ->
immutable QuCN{B<:AbstractBasis,T,A} <: QuantumPropagator{B,T,A}
    hamiltonian::QuMatrix{B,T,A}
    init_state::QuVector
    tlist::Array
    QuCN(hamiltonian, init_state, tlist) = new(hamiltonian, init_state, tlist)
end

@doc """
Crank Nicolson State
`psi`     : State of the system
`t_index` : Array index of time in `tlist` of the current state of the system
""" ->
immutable QuCNState
    psi::QuVector
    t_index::Int
end
QuCN{B<:AbstractBasis,T,A}(hamiltonian::QuMatrix{B,T,A}, init_state::QuVector, tlist::Array)=QuCN{B,T,A}(hamiltonian, init_state, tlist)

@doc """
Input Parameters : Crank Nicolson Problem

Returns the starting, i.e., the initial state of the system
""" ->
function Base.start(prob::QuCN)
    return QuCNState(prob.init_state, 1)
end

@doc """
Input Parameters : Crank Nicolson Problem and Crank Nicolson State

Returns the next state using Euler's method.
""" ->
function Base.next(prob::QuCN, state::QuCNState)
    op = prob.hamiltonian
    t_list = prob.tlist
    current_state = state.psi
    current_t = state.t_index
    if t_list[current_t] != t_list[end]
        next_t = t_list[current_t+1]
        uni = eye(op)-im*op*(next_t-t_list[current_t])/2
        next_state =  \(uni', uni*current_state)
        return (next_state, next_t), QuCNState(next_state, findfirst(t_list, next_t))
    else
        return QuCNState(current_state, findfirst(t_list,current_t))
    end
end

@doc """
Input Parameters : Crank Nicolson Problem and Crank Nicolson State

Returns true if the current state is final state, else false
""" ->
function Base.done(prob::QuCN, state::QuCNState)
     if prob.tlist[state.t_index] == prob.tlist[end]
       return true
     else
       return false
     end
end

@doc """
Input Parameters : Hamiltonian of the system, Initial State of the system, Time steps

Returns the collection of output states at every time step using Crank Nicolson Method.
""" ->
function solver_CN(hamiltonian::QuMatrix, init_state::QuVector, tlist::Array)
    final_state = initialize(hamiltonian, init_state, tlist)
    #for i=2:length(tlist)
    #    uni = QuArray(eye(size(hamiltonian,1)))-(im*hamiltonian*(tlist[i]-tlist[i-1]))/2
    #    push!(final_state, \(uni',uni*final_state[i-1]))
    #end
    id = eye(hamiltonian)
    next_step = start(tlist)
    while !done(tlist, next_step)
        step, next_step = next(tlist, next_step)
        if next_step == length(tlist)+1
          break
        else
          uni = id-(im*hamiltonian*(tlist[next_step]-step))/2
          push!(final_state, \(uni',uni*final_state[next_step-1]))
        end
    end
    return final_state
end

export QuEuler,
      QuEulerState,
      solver_euler,
      solver_CN
      start,
      next,
      done