from sympy import *
  from sympy.simplify.fu import *
  
  init_printing()
  
  #constants
  E0, w_laser, W_aom, t = symbols('E0, omega_laser, Omega_aom, t', real=True, imaginary=False)
  dphi, dW = symbols('Delta_phi, Delta_Omega', real=True, imaginary=False)
  c, dx = symbols('c, Delta_x', real=True, imaginary=False)
  
  phi_u = dW*t + dphi
  phi = W_aom*t + phi_u
  
  #laser
  E_laser = E0*exp(I*(w_laser*t))
  
  #reflected wave and dephased by doppler and noise and double shifted
  E_aom_ref = E0*exp(I*( w_laser*t + 2*dx*(w_laser+phi.diff(t))/c + 2*phi ))
  
  #reflected wave on the reference arm of michelson
  E_ref = E_laser
  
  #optical mixer
  E_mix = sqrt(2)/2*E_ref + sqrt(2)/2*E_aom_ref
  
  #intensity of mixed waves
  I_mix = abs(E_mix)**2
  I_mix = expand(TR8(expand(expand_complex(I_mix))))
  
  #Q demodulation
  kphi = symbols('k_phi', real=True, imaginary=False)
  #eps = 2 * kphi * I_mix * cos(2*phi.diff(t)*t-pi/2)
  eps = 2 * kphi * I_mix * cos(2*W_aom*t-pi/2)
  eps = expand(TR8(TR7(expand(eps))))