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))