arrayintg_toy.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
##
# arrayintg_toy.py: Toy program to illustrate problems integrating an array with SciPy.
##
# © 2010 Christopher E. Granade (cgranade@gmail.com).
# Licensed under the GPL version 2.
##
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
##

## IMPORTS #####################################################################

import math as m
import numpy as np
import scipy.linalg.matfuncs as mf

import scipy.integrate

## CONSTANTS ###################################################################

## OTHER ##
I = 1j

## COMMON MATRICES ##
Id = np.array([[1,0],[0,1]])
sX = np.array([[0,1],[1,0]])
sY = np.array([[0,-I],[I,0]])
sZ = np.array([[1,0],[0,-1]])

Ep = np.array([[1,0],[0,0]])
Em = np.array([[0,0],[0,1]])
sp = np.array([[0,1],[0,0]])
sm = np.array([[0,0],[1,0]])

## FUNCTIONS ###################################################################

def propagator(H, t):
    return mf.expm(-I * t * H)

def interaction_frame(H0, Hint, t):
    return reduce(np.dot,[
        propagator(Hint, t),
        H0,
        propagator(Hint, -t)
    ])

def H0(omega):
    return omega * np.kron(sZ, sZ)

## MAIN ########################################################################

if __name__ == "__main__":

    wz =  100.0
    wt = 1000.0

    H_cycle = (m.pi / wt) * scipy.integrate.romberg(
       lambda t: interaction_frame(
            H0(wz),
            wt * (np.kron(sX, Id) + np.kron(Id, sX)),
            t
        ),
        0, m.pi / wt
    )

    print H_cycle