def rmpHomozygousWithTheta(theta, pu):
  """
  Notation from Baldings and Nichols (1994)
  taken from Buckleton FSI Gen (23) 2016; 91-100
  https://doi.org/10.1016/j.fsigen.2016.03.004
  """
  # 2theta + (1-theta)*pu
  numerator1 = (theta+theta) + (1-theta)*pu
  # 3theta + (1-theta)*pu == theta + numerator1
  numerator2 = theta + numerator1
  
  # 1 + theta TIMES 1 + 2theta
  denom = (1+theta) * (1 + theta + theta)

  return (numerator1*numerator2)/denom