def Read_in_data(scifile,varfile,SURVEY='PS1',vb=False):

   hdulist = pyfits.open(scifile)
   sci = hdulist[0].data
   hdr = hdulist[0].header
   hdulist.close()
   NX,NY = sci.shape

   if (varfile is not None): 
       hdulist = pyfits.open(varfile)
       var = hdulist[0].data
       hdulist.close()
   else:
   # Make a var image from the sci image...
       background = np.median(sci)
       diffimage = sci - background
   # Get the flux-to-count conversion factor from header, at least in SDSS:
       try:
          tmpsurvey = hdr['ORIGIN']
          if (tmpsurvey == 'SDSS'):
             tempmtoc = hdr['NMGY']
          else:
             tempmtoc = 1.
       except:
          tempmtoc = 1.