docs: def define_laser_variables(Nl, real_amplitudes=False, variables=None):
      r"""Return the amplitudes and frequencies of Nl fields.
  
      >>> E0, omega_laser = define_laser_variables(2)
      >>> E0, omega_laser
      ([E_0^1, E_0^2], [varpi_1, varpi_2])
  
      The amplitudes are complex by default:
      >>> conjugate(E0[0])
      conjugate(E_0^1)
  
      But they can optionally be made real:
      >>> E0, omega_laser = define_laser_variables(2, real_amplitudes=True)
      >>> conjugate(E0[0])
      E_0^1
  
      They can also be made explicit functions of given variables:
      >>> from sympy import symbols
      >>> t, z = symbols("t, z", real=True)
      >>> E0, omega_laser = define_laser_variables(2, variables=[t, z])
      >>> E0
      [E_0^1(t, z), E_0^2(t, z)]
  
      """
      if variables is None:
          E0 = [Symbol(r"E_0^"+str(l+1), real=real_amplitudes)
                for l in range(Nl)]
      else:
          E0 = [Function(r"E_0^"+str(l+1), real=real_amplitudes)(*variables)
                for l in range(Nl)]
  
      omega_laser = [Symbol(r"varpi_"+str(l+1), positive=True)
                     for l in range(Nl)]
      return E0, omega_laser