docs: def _initialize_background(self):
          """Set up background state (zonal flow and PV gradients)."""
  
          # Background zonal flow (m/s):
          self.H = self.Hi.sum()
          self.set_U1U2(self.U1, self.U2)
          self.U = self.U1 - self.U2
  
          # the F parameters
          self.F1 = self.rd**-2 / (1.+self.delta)
          self.F2 = self.delta*self.F1
  
          # the meridional PV gradients in each layer
          # need to calculate actual PV gradient
          self.Qy1 = self.beta + self.F1*(self.U1 - self.U2) + np.gradient(np.gradient(self.U1,self.dy),self.dy)
          self.Qy2 = self.beta - self.F2*(self.U1 - self.U2) + np.gradient(np.gradient(self.U2,self.dy),self.dy)
          self.Qy = np.array([self.Qy1, self.Qy2])
          # complex versions, multiplied by k, speeds up computations to precompute
          self.ikQy1 = self.Qy1[:,np.newaxis] * 1j * self.k
          self.ikQy2 = self.Qy2[:,np.newaxis] * 1j * self.k
  
          # vector version
          self.ikQy = np.vstack([self.ikQy1[np.newaxis,...],
                                 self.ikQy2[np.newaxis,...]])
          self.ilQx = 0.
  
          # layer spacing
          self.del1 = self.delta/(self.delta+1.)
          self.del2 = (self.delta+1.)**-1