from sympy import *
from sympy.simplify.fu import *

init_printing()

#constants
E0, E1, w, W, t, dphi, kphi = symbols('E0, E1, omega, Omega, t, Delta_phi, k_phi', real=True, imaginary=False)
dx, du, c = symbols('Delta_x, Delta_u, c', real=True, imaginary=False)

#reference wave
Eb0 = E0*exp(I*(w*t))

#double Omega shifted and dphi dephased wave
Eb1 = E1*exp(I*((w+2*W)*t -2*w*(dx+2*du)/c))

#optical mixer
Emix = sqrt(2)/2 * Eb0 + sqrt(2)/2 * Eb1

#intensity of mixed waves
Imix = simplify(expand_complex(abs(Emix)**2))

#demodulation Q at 2 * Omega
Vmix = 2 * kphi * Imix * cos(2*W*t-pi/2)
Vmix = expand(TR8(TR7(expand(Vmix))))
print(pretty(Vmix))