docs: def liquid_jet_pump_ancillary(rhop, rhos, Kp, Ks, d_nozzle=None, d_mixing=None, 
                                Qp=None, Qs=None, P1=None, P2=None):
      r'''Calculates the remaining variable in a liquid jet pump when solving for
      one if the inlet variables only and the rest of them are known. The
      equation comes from conservation of energy and momentum in the mixing
      chamber.
      
      The variable to be solved for must be one of `d_nozzle`, `d_mixing`,
      `Qp`, `Qs`, `P1`, or `P2`.
      
      .. math::
          P_1 - P_2 = \frac{1}{2}\rho_pV_n^2(1+K_p)
          - \frac{1}{2}\rho_s V_3^2(1+K_s)
          
      Rearrange to express V3 in terms of Vn, and using the density ratio `C`,  
      the expression becomes:
          
      .. math::
          P_1 - P_2 = \frac{1}{2}\rho_p V_n^2\left[(1+K_p) - C(1+K_s) 
          \left(\frac{MR}{1-R}\right)^2\right]
  
      Using the primary nozzle area and flow rate:
          
      .. math::
          P_1 - P_2 = \frac{1}{2}\rho_p \left(\frac{Q_p}{A_n}\right)^2
          \left[(1+K_p) - C(1+K_s) \left(\frac{MR}{1-R}\right)^2\right]
  
      For `P`, `P2`, `Qs`, and `Qp`, the equation can be rearranged explicitly 
      for them. For `d_mixing` and `d_nozzle`, a bounded solver is used searching
      between 1E-9 m and 20 times the other diameter which was specified.
  
      Parameters
      ----------
      rhop : float
          The density of the primary (motive) fluid, [kg/m^3]
      rhos : float
          The density of the secondary fluid (drawn from the vacuum chamber),
          [kg/m^3]
      Kp : float
          The primary nozzle loss coefficient, [-]
      Ks : float
          The secondary inlet loss coefficient, [-]
      d_nozzle : float, optional
          The inside diameter of the primary fluid's nozle, [m]
      d_mixing : float, optional
          The diameter of the mixing chamber, [m]
      Qp : float, optional
          The volumetric flow rate of the primary fluid, [m^3/s]
      Qs : float, optional
          The volumetric flow rate of the secondary fluid, [m^3/s]
      P1 : float, optional
          The pressure of the primary fluid entering its nozzle, [Pa]
      P2 : float, optional
          The pressure of the secondary fluid at the entry of the ejector, [Pa]
          
      Returns
      -------
      solution : float
          The parameter not specified (one of `d_nozzle`, `d_mixing`,
          `Qp`, `Qs`, `P1`, or `P2`), (units of `m`, `m`, `m^3/s`, `m^3/s`,
          `Pa`, or `Pa` respectively)
  
      Notes
      -----
      The following SymPy code was used to obtain the analytical formulas (
      they are not shown here due to their length):
      
      >>> from sympy import *
      >>> A_nozzle, A_mixing, Qs, Qp, P1, P2, rhos, rhop, Ks, Kp = symbols('A_nozzle, A_mixing, Qs, Qp, P1, P2, rhos, rhop, Ks, Kp')
      >>> R = A_nozzle/A_mixing
      >>> M = Qs/Qp
      >>> C = rhos/rhop
      >>> rhs = rhop/2*(Qp/A_nozzle)**2*((1+Kp) - C*(1 + Ks)*((M*R)/(1-R))**2 )
      >>> new = Eq(P1 - P2,  rhs)
      >>> #solve(new, Qp)
      >>> #solve(new, Qs)
      >>> #solve(new, P1)
      >>> #solve(new, P2)
  
      Examples
      --------
      Calculating primary fluid nozzle inlet pressure P1:
      
      >>> liquid_jet_pump_ancillary(rhop=998., rhos=1098., Ks=0.11, Kp=.04, 
      ... P2=133600, Qp=0.01, Qs=0.01, d_mixing=0.045, d_nozzle=0.02238)
      426434.60314398084
  
      References
      ----------
      .. [1] Ejectors and Jet Pumps. Design and Performance for Incompressible 
         Liquid Flow. 85032. ESDU International PLC, 1985.
      '''
      unknowns = sum(i is None for i in (d_nozzle, d_mixing, Qs, Qp, P1, P2))
      if unknowns > 1:
          raise Exception('Too many unknowns')
      elif unknowns < 1:
          raise Exception('Overspecified')
      C = rhos/rhop
      
      if Qp is not None and Qs is not None:
          M = Qs/Qp
      if d_nozzle is not None:
          A_nozzle = pi/4*d_nozzle*d_nozzle
          if d_mixing is not None:
              A_mixing = pi/4*d_mixing*d_mixing
              R = A_nozzle/A_mixing
              
      if P1 is None:        
          return rhop/2*(Qp/A_nozzle)**2*((1+Kp) - C*(1 + Ks)*((M*R)/(1-R))**2 ) + P2
      elif P2 is None:
          return -rhop/2*(Qp/A_nozzle)**2*((1+Kp) - C*(1 + Ks)*((M*R)/(1-R))**2 ) + P1
      elif Qs is None:
          try:
              return ((-2*A_nozzle**2*P1 + 2*A_nozzle**2*P2 + Kp*Qp**2*rhop + Qp**2*rhop)/(C*rhop*(Ks + 1)))**0.5*(A_mixing - A_nozzle)/A_nozzle
          except ValueError:
              return -1j
      elif Qp is None:
          return A_nozzle*((2*A_mixing**2*P1 - 2*A_mixing**2*P2 - 4*A_mixing*A_nozzle*P1 + 4*A_mixing*A_nozzle*P2 + 2*A_nozzle**2*P1 - 2*A_nozzle**2*P2 + C*Ks*Qs**2*rhop + C*Qs**2*rhop)/(rhop*(Kp + 1)))**0.5/(A_mixing - A_nozzle)
      elif d_nozzle is None:
          def err(d_nozzle):
              return P1 - liquid_jet_pump_ancillary(rhop=rhop, rhos=rhos, Kp=Kp, Ks=Ks, d_nozzle=d_nozzle, d_mixing=d_mixing, Qp=Qp, Qs=Qs, 
                                P1=None, P2=P2)
          return brenth(err, 1E-9, d_mixing*20)
      elif d_mixing is None:
          def err(d_mixing):
              return P1 - liquid_jet_pump_ancillary(rhop=rhop, rhos=rhos, Kp=Kp, Ks=Ks, d_nozzle=d_nozzle, d_mixing=d_mixing, Qp=Qp, Qs=Qs, 
                                P1=None, P2=P2)
          try:
              return brenth(err, 1E-9, d_nozzle*20)
          except:
              return newton(err, d_nozzle*2)