from sympy import *
init_printing()

# Note: This notebook requires and experimental branch of sympy available at: https://github.com/jrjohansson/sympy
from sympy.physics.quantum import *

theta = Symbol("theta", real=True)
x = Symbol("x", real=True)

dL, v, wp, wm = symbols("{\delta\!L}, v, omega_-, omega_+")

# input field operators
am = BosonOperator("a_-")
ap = BosonOperator("a_+")

# output field operators
bm = BosonOperator("b_-")
bp = BosonOperator("b_+")

# relation between input field operators and output field operators
bm_sub = -am - 1j * x * Dagger(ap)
bp_sub = -ap - 1j * x * Dagger(am)

bmd_sub = Dagger(bm_sub)
bpd_sub = Dagger(bp_sub)

f = exp(1j*theta) * bm + exp(-1j*theta) * Dagger(bm) + 1j * (exp(1j*theta) * bp - exp(-1j*theta) * Dagger(bp))

f

fdf = Dagger(f) * f

fdf.expand()

fdf_normal = normal_order(fdf.expand())

fdf_normal

fdf_1 = fdf_normal.expand().subs({bm: bm_sub, Dagger(bm): bmd_sub, bp: bp_sub, Dagger(bp): bpd_sub}).expand()

fdf_1

fdf_2 = normal_ordered_form(fdf_1.expand(), independent=True)

fdf_2

# assume thermal input: first remove nonthermal operator combinations
fdf_3 = fdf_2.subs({am**2: 0, ap**2: 0,
                    ap*am: 0, am*ap: 0, 
                    Dagger(ap) * am: 0, Dagger(am) * ap: 0, 
                    Dagger(ap) * Dagger(am): 0, Dagger(am) * Dagger(ap): 0, 
                    Dagger(ap)**2: 0, 
                    Dagger(am)**2: 0})

fdf_3

n_th_m, n_th_p = symbols("{n^-_{th}}, {n^+_{th}}", positive=True)

# thermal operator combinations 
fdf_4 = fdf_3.subs({Dagger(am) * am: n_th_m, Dagger(ap) * ap: n_th_p})

fdf_4

fdf_5 = fdf_4.collect(x).replace(exp(2j*theta), cos(2*theta) + 1j * sin(2*theta)).\
                         replace(exp(-2j*theta), cos(2*theta) - 1j * sin(2*theta))
    
fdf_6 = fdf_5.simplify().collect(cos(2*theta)).collect(x).simplify()

fdf_6

fdf_7 = fdf_6.subs(x**2, 0)

fdf_7

fdf_8 = fdf_7.subs(x**2, 0).subs(x, dL/v * sqrt(wp * wm))

fdf_8

# switch back to using x as small paramter, and set theta = 0
fdf_9 = fdf_8.subs(dL/v * sqrt(wp*wm), x).subs(theta, 0)

fdf_9

fdf_neg = solve(fdf_9, x)[0]

fdf_neg

#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%reload_ext version_information

%version_information sympy