code: def covariance_matrix(self, data, background_rms=1, exposure_map=1, noise_map=None, verbose=False):
          
          if noise_map is not None:
              return noise_map**2
          if isinstance(exposure_map, int) or isinstance(exposure_map, float):
              if exposure_map <= 0:
                  exposure_map = 1
          else:
              mean_exp_time = np.mean(exposure_map)
              exposure_map[exposure_map < mean_exp_time / 10] = mean_exp_time / 10
          if verbose:
              if background_rms * np.max(exposure_map) < 1:
                  print("WARNING! sigma_b*f %s < 1 count may introduce unstable error estimates" % (background_rms * np.max(exposure_map)))
          d_pos = np.zeros_like(data)
          #threshold = 1.5*sigma_b
          d_pos[data >= 0] = data[data >= 0]
          #d_pos[d < threshold] = 0
          sigma = d_pos / exposure_map + background_rms ** 2
          return sigma
docs: returns a diagonal matrix for the covariance estimation which describes the error
  
          Notes:
  
          - the exposure map must be positive definite. Values that deviate too much from the mean exposure time will be
              given a lower limit to not under-predict the Poisson component of the noise.
  
          - the data must be positive semi-definite for the Poisson noise estimate.
              Values < 0 (Possible after mean subtraction) will not have a Poisson component in their noise estimate.
  
  
          :param data: data array, eg in units of photons/second
          :param background_rms: background noise rms, eg. in units (photons/second)^2
          :param exposure_map: exposure time per pixel, e.g. in units of seconds
          :return: len(d) x len(d) matrix that give the error of background and Poisson components; (photons/second)^2