In [1]:

    

import numpy as np
import math
import matplotlib.pyplot as plt

import scipy
import scipy.optimize 

i = np.complex(0, 1)


    

In [63]:

    

""" Written-down solution of Schinger's spin-half Hamiltonian. """

def Eval1 (t, w0):
    return -w0/2.

def Eval2 (t, w0):
    return w0/2.

def Evect1 (t, w, w0, theta):
    return np.array([ np.exp(-i * w * t/2.) * np.sin(theta/2.),
                     -np.exp(i * w * t/2.)*np.cos(theta/2.)])

def Evect2 (t, w, w0, theta):
    return np.array([ np.exp(-i * w * t/2.) * np.cos(theta/2.),
                     np.exp(i * w * t/2.)*np.sin(theta/2.)])

def wtilde (w, w0, theta):
    return np.sqrt(w0**2 + w**2 - 2*w0*w*np.cos(theta))

def c1 (t, w, w0, theta):
    wt = wtilde(w, w0, theta)
    return np.cos(wt/2.) + i*np.sin(wt/2.)*(w0-w*np.cos(theta))/wt

def c2 (t, w, w0, theta):
    wt = wtilde(w, w0, theta)
    return i*(w/wt)*np.sin(theta)*np.sin(wt*t/2.)

def beta1 (t, w, w0, theta):
    return w0*t/2. - (w*t/2.)*np.cos(theta)


# Terms that contribute to (9).
def Q2 (t, w, w0, theta, beta):
    return np.exp(i*beta)*(w/(2.*w0))*np.sin(theta)

def R2 (t, w, w0, theta, beta):
    return c2(t, w, w0, theta) - Q2(t, w, w0, theta, beta)


    

In [40]:

    

t = np.linspace(0, 50, 700)


    

In [20]:

    

def max_c2 (w, w0, theta):
    """
        Helper function to determine the maximum function value of c2 w.r.t t, numerically.
    """
    max_t = scipy.optimize.fmin(lambda t: -np.abs(c2(t, w, w0, theta)), 0, disp=False)
    return np.abs(c2(max_t[0], w, w0, theta))


    

In [64]:

    

def plotInto (t, p, w, w0, theta, xlim, ylim):
    c2_max = max_c2(w, w0, theta)
    p.plot(t, np.abs(c2(t/w0, w, w0, theta)), "b", label=r"$|c_2|$")
    p.plot(t, np.abs(Q2(t/w0, w, w0, theta, beta1(t/w0, w, w0, theta))), "r", label=r"$|Q_2$|")
    p.plot(t, np.abs(R2(t/w0, w, w0, theta, beta1(t/w0, w, w0, theta))), "g", label=r"$|R_2|$")
    p.plot(t, [c2_max for x in xrange(len(t))], color="black", ls="--",
           label="max $c_2 = %0.2f$" % np.round(c2_max, 2))
    p.set_xlabel(r'Time in units of $1/w_0$.', fontsize=18)
    p.set_title(r"$w/w_0=%0.1f,\ \theta=%0.2f\times \pi/2$" % (w/float(w0), theta*(2./np.pi)))
    p.set_xlim(xlim)
    p.set_ylim(ylim)
    p.legend(loc=2)
    

fig, ax = plt.subplots(1, 2, figsize=(20, 6))

plotInto(t, ax[0], 1, 10., 0.1*pion2 , [0, 50], [0, 0.12])
plotInto(t, ax[1], 10., 1., 0.1*pion2, [0, 5], [0, 6])

fig, ax = plt.subplots(1, 2, figsize=(20, 6))
plotInto(t, ax[0], 1., 10., 0.5*pion2, [0, 50], [0, 0.12])
plotInto(t, ax[1], 10., 1., 0.5*pion2, [0, 5], [0, 6])

fig, ax = plt.subplots(1, 2, figsize=(20, 6))
plotInto(t, ax[0], 1., 10., 1*pion2, [0, 50], [0, 0.12])
plotInto(t, ax[1], 10., 1., 1*pion2, [0, 5], [0, 6])