def Initial_PSF(FWHM,double=False):

   # NB. FWHM of PSF is given in pixels.

   if not double:
      # Single Gaussian default:
      w = np.array([1.0])                      # amplitude at peak
      mu = np.array([[0.0,0.0]])               # centroid position in pixels 
      var = (FWHM/2.35)**2.0
      cov = np.array([[[var,0.0],[0.0,var]]])  # pixels^2, covariance matrix

   else:
      # Double Gaussian alternative:
      w = np.array([1.0,1.0])                                        
      mu = np.array([[0.0,0.0],[0.0,0.0]])                           
      var = (FWHM/2.35)**2.0
      cov = np.array([[[1.0,0.0],[0.0,1.0]],[[var,0.0],[0.0,var]]])