code: def restore(self):
          
          clean_beam, beam_params = beam_fit(self.psf_data, self.cdelt1, self.cdelt2)
  
          if np.all(np.array(self.psf_data_shape)==2*np.array(self.dirty_data_shape)):
              self.restored = np.fft.fftshift(np.fft.irfft2(np.fft.rfft2(conv.pad_array(self.model))*np.fft.rfft2(clean_beam)))
              self.restored = self.restored[self.dirty_data_shape[0]/2:-self.dirty_data_shape[0]/2,
                              self.dirty_data_shape[1]/2:-self.dirty_data_shape[1]/2]
          else:
              self.restored = np.fft.fftshift(np.fft.irfft2(np.fft.rfft2(self.model)*np.fft.rfft2(clean_beam)))
          self.restored += self.residual
          self.restored = self.restored.astype(np.float32)
  
          return beam_params
docs: This method constructs the restoring beam and then adds the convolution to the residual.