diff --git a/src/simmer/analyze_image.py b/src/simmer/analyze_image.py
index 9a03c86..ae83919 100644
--- a/src/simmer/analyze_image.py
+++ b/src/simmer/analyze_image.py
@@ -27,6 +27,9 @@
 
 import simmer.contrast as sim_con_curve
 
+import logging
+logger = logging.getLogger('simmer')
+
 #Test file
 startname = '/Users/courtney/Documents/data/shaneAO/data-'
 datestr = '2019-07-19'
@@ -39,6 +42,8 @@
 outdir = startname + datestr + midname + datestr + '/'+starname+'/'+filt+'/'
 
 def analyze(filename=filename, maxiter = 10, postol=1, fwhmtol = 0.5, inst = 'ShARCS', outdir='', verbose=True):
+    if verbose:
+        logger.setLevel(logging.DEBUG)
 
     #set defaults
     if inst == 'PHARO':
@@ -58,8 +63,8 @@ def analyze(filename=filename, maxiter = 10, postol=1, fwhmtol = 0.5, inst = 'Sh
 
     #Determine FWHM using that center
     fwhm = find_FWHM(im, [xcen,ycen])
-    if verbose == True:
-        print('Estimated FWHM: ', fwhm)
+    logger.debug('Estimated FWHM: ', fwhm)
+
 
     #Iterate until coordinates and FWHM agree to within tolerance
     #or until maximum number of iterations is reached
@@ -67,13 +72,15 @@ def analyze(filename=filename, maxiter = 10, postol=1, fwhmtol = 0.5, inst = 'Sh
     fwhmdiff = 2
     niter = 1
     while np.logical_and(niter < maxiter, np.logical_or(posdiff > postol, fwhmdiff > fwhmtol)):
-        if verbose == True:
-            print('Beginning iteration ', niter)
+
+        logger.debug('Beginning iteration ', niter)
+
 
         #Find sources again
         updated_sources = find_sources(im, fwhm=fwhm)
-        if verbose == True:
-            print('Updated sources')
+        logger.debug('Updated sources')
+
+        if verbose:
             print(updated_sources)
 
         #Find brightest peak again using updated list of stars
@@ -81,15 +88,16 @@ def analyze(filename=filename, maxiter = 10, postol=1, fwhmtol = 0.5, inst = 'Sh
 
         #Determine FWHM using that center
         updated_fwhm = find_FWHM(im, [updated_xcen,updated_ycen])
-        if verbose == True:
-            print('Estimated FWHM: ', updated_fwhm)
+        logger.debug('Estimated FWHM: ', updated_fwhm)
+
 
         #Compute differences
         posdiff = np.sqrt((updated_xcen - xcen)**2. + (updated_ycen - ycen)**2.)
         fwhmdiff = np.sqrt((updated_fwhm - fwhm)**2.)
-        if verbose == True:
-            print('Current posdiff: ', posdiff)
-            print('Current fwhmdiff: ', fwhmdiff)
+
+        logger.debug('posdiff: ', posdiff)
+        logger.debug('fwhmdiff: ', fwhmdiff)
+
 
         #Update reference values
         xcen = updated_xcen
@@ -117,21 +125,23 @@ def find_sources(im, sigma=5, fwhm=5, tscale=10, verbose=False, plot=True, **kwa
     """
     Determines sources in an image. Based on astropy tutorial here: https://photutils.readthedocs.io/en/stable/detection.html
     """
+    if verbose:
+        logger.setLevel(logging.DEBUG)
     mean, median, std = sigma_clipped_stats(im, sigma=sigma)
-    if verbose == True:
-        print((mean, median, std))
+    logger.debug('mean, median, std: ', mean, median, std)
+
 
     sources = None
     while np.logical_and(type(sources) == type(None), fwhm < 200):
-        if verbose == True:
-            print('trying FWHM: ', fwhm)
+        logger.debug('trying FWHM: ', fwhm)
+
         #Detect stars >threshold above the background. Needed to adjust FWHM and threshold
         daofind = DAOStarFinder(fwhm=fwhm, threshold=tscale*std, exclude_border=True, **kwargs)
         sources = daofind(im - median)
     #    for col in sources.colnames:
     #        sources[col].info.format = '%.8g'  # for consistent table output
-        if verbose == True:
-            print(sources)
+        logger.debug('sources: ', sources)
+
         fwhm += 1
 
     if plot:
@@ -145,7 +155,6 @@ def find_sources(im, sigma=5, fwhm=5, tscale=10, verbose=False, plot=True, **kwa
 
     #Convert sources to a dataframe
     df = sources.to_pandas()
-   # print(len(sources))
     return df
 
 def find_FWHM(image, center, min_fwhm = 2, verbose=False):
@@ -208,9 +217,12 @@ def aperture_photometry(im, df, fwhm):
 
 def find_center(df, verbose=False):
     '''Returns coordinates of star with highest peak'''
+    if verbose:
+        logger.setLevel(logging.DEBUG)
+
     ww = np.where(df['peak'] == np.max(df['peak']))
     xcen = int(np.round(float(df.iloc[ww]['xcentroid'])))
     ycen = int(np.round(float(df.iloc[ww]['ycentroid'])))
-    if verbose == True:
-        print('Closest pixels to center: ', xcen, ycen)
+    logger.debug('Closest pixels to center: ', xcen, ycen)
+
     return xcen, ycen
diff --git a/src/simmer/check_logsheet.py b/src/simmer/check_logsheet.py
index 2558dcb..f282e65 100644
--- a/src/simmer/check_logsheet.py
+++ b/src/simmer/check_logsheet.py
@@ -11,6 +11,8 @@
 import pandas as pd
 
 from . import add_dark_exp as ad
+import logging
+logger = logging.getLogger('simmer')
 
 
 def check_logsheet(inst, log_name, tab=None, add_dark_times=False):
@@ -57,58 +59,60 @@ def check_tab(inst, add_dark_times, tab=None):
             if not np.isin(col, frame_cols):
                 missing.append(col)
         if len(missing) != 0:
-            print(f"Missing columns for {missing}.")
+            logger.error(f"Missing columns for {missing}.")
             failed += 1
 
         objects = log_frame["Object"].dropna().values
         exptimes = log_frame["ExpTime"].dropna().values
         if len(exptimes) != len(objects):
-            print("Missing an exposure time.")
+            logger.error("Missing an exposure time.")
             failed += 1
 
         filters = log_frame["Filter"].dropna().values
         if len(filters) != len(log_frame[log_frame["Object"] != "dark"]):
-            print("Missing a filter.")
+            logger.error("Missing a filter.")
             failed += 1
 
         starts = log_frame["Start"].dropna().values
         if len(starts) != len(objects):
-            print("Missing a start exposure.")
+            logger.error("Missing a start exposure.")
             failed += 1
 
         ends = log_frame["End"].dropna().values
         if len(ends) != len(objects):
-            print("Missing an end exposure.")
+            logger.error("Missing an end exposure.")
             failed += 1
 
         coadds = log_frame["Coadds"].dropna().values
         if len(coadds) != len(objects):
-            print("Missing a coadd.")
+            logger.error("Missing a coadd.")
+
             failed += 1
 
         if not np.all(exptimes > 0):
-            print("There are negative or 0 exposure times.")
+            logger.error("There are negative or 0 exposure times.")
+
             failed += 1
         try:
             inter = ends - starts
             if not np.all(inter >= 0):
-                print("Check the start and end exposures.")
+                logger.error("There are negative exposure times.")
+
                 failed += 1
         except ValueError:
-            print("Check the start and end exposures.")
+            logger.error("Check the start and end exposures.")
+
             failed += 1
 
         exposes = log_frame["Expose"].dropna().values
         try:
             if not np.all(exposes == inter + 1):
-                print(
-                    "Incorrect number of exposures for start and end exposure."
-                )
+                logger.error('Incorrect number of exposures for start and end exposure.')
                 failed += 1
         except UnboundLocalError:
-            print("Incorrect number of exposures for start and end exposure.")
+            logger.error('Incorrect number of exposures for start and end exposure.')
             failed += 1
-        print(f"{9-failed}/9 logsheet checks passed.")
+        logger.info(f"{9-failed}/9 logsheet checks passed.")
         return failed
 
     failed = 0
diff --git a/src/simmer/drivers.py b/src/simmer/drivers.py
index d18390e..e660ffd 100644
--- a/src/simmer/drivers.py
+++ b/src/simmer/drivers.py
@@ -15,6 +15,10 @@
 from . import sky
 from . import summarize as summarize
 
+import logging
+
+logger = logging.getLogger('simmer')
+
 def all_driver(
 
     inst, config_file, raw_dir, reddir, sep_skies = False, plotting_yml=None, searchsize=10, just_images=False, verbose=True
@@ -77,8 +81,8 @@ def all_driver(
             leave=True,
         )
     ):
-        print('searchsize: ', searchsize)
-        print('s_dir: ', s_dir)
+        logger.info(f"Running registration for {s_dir}")
+        logger.info(f"searchsize: {searchsize}")
         image.create_im(s_dir, searchsize, method=methods[i], verbose=verbose)
 
     #make summary plot showing reduced images of all stars observed
diff --git a/src/simmer/image.py b/src/simmer/image.py
index 58e745a..2e1314a 100644
--- a/src/simmer/image.py
+++ b/src/simmer/image.py
@@ -17,6 +17,9 @@
 from . import utils as u
 from . import contrast as contrast
 
+import logging
+logger = logging.getLogger('simmer')
+
 
 class FlatOpeningError(ValueError):
     pass
@@ -254,12 +257,10 @@ def create_im(s_dir, ssize1, plotting_yml=None, fdirs=None, method="quick_look",
         #Only register star images, not sky images
         dirparts = sf_dir.split('/')
         if 'sky' in dirparts[len(dirparts)-3]:
-            if verbose == True:
-                print('this is a sky directory: ', sf_dir)
+            logger.debug('this is a sky directory: ', sf_dir)
             continue
 
-        if verbose == True:
-            print('working on sf_dir ', sf_dir)
+        logger.debug('working on sf_dir ', sf_dir)
 
         files = glob(
             sf_dir + f"sh*.fits"
@@ -300,7 +301,7 @@ def create_im(s_dir, ssize1, plotting_yml=None, fdirs=None, method="quick_look",
                 image[image < 0.0] = 0.0
                 image_centered = reg.register_bruteforce(image)
                 if len(image_centered) == 0:
-                    print("Resorting to saturated mode.")
+                    logger.info("Resorting to saturated mode.")
                     image_centered, rot, newshifts1 = reg.register_saturated(
                         image, ssize1, newshifts1
                     )
@@ -329,9 +330,9 @@ def create_im(s_dir, ssize1, plotting_yml=None, fdirs=None, method="quick_look",
         aend = astart+cutsize
         bend = bstart+cutsize
         if np.logical_or(aend > final_im.shape[0],bend > final_im.shape[1]):
-            print('ERROR: Requested cutout is too large. Using full image instead.')
-            print('Current image dimensions: ', final_im.shape)
-            print('Desired cuts: ', astart, aend, bstart, bend)
+            logger.error('ERROR: Requested cutout is too large. Using full image instead.')
+            logger.info('Current image dimensions: ', final_im.shape)
+            logger.info('Desired cuts: ', astart, aend, bstart, bend)
         else:
             final_im = final_im[astart:astart+cutsize,bstart:bstart+cutsize] #extract central cutsize x cutsize pixel region from larger image
 
diff --git a/src/simmer/make_triceratops_contrasts.py b/src/simmer/make_triceratops_contrasts.py
index 65ae2d5..c4816ee 100644
--- a/src/simmer/make_triceratops_contrasts.py
+++ b/src/simmer/make_triceratops_contrasts.py
@@ -12,72 +12,82 @@
 import pandas as pd
 import os as os
 import glob as glob
+import logging
 
-#Defaults
-verbose = False
+logger = logging.getLogger('simmer')
 
+if __name__ == '__main__':
 
-#Files and directories
-tridir = '/Users/courtney/Documents/data/toi_paper_data/contrast_curves_for_triceratops/'
+    #Defaults
+    verbose = False
 
-#Get list of final images
-namestr = '/Users/courtney/Documents/data/shaneAO/*/reduced*/*/*/contrast_curve.csv'
-flist = glob.glob(namestr)
+    if verbose:
+        logger.setLevel(logging.DEBUG)
 
-print('Files: ', len(flist))
 
-#Loop through these files!
-counter = 0
-for ff in np.arange(len(flist)):
+    #Files and directories
+    tridir = '/Users/courtney/Documents/data/toi_paper_data/contrast_curves_for_triceratops/'
 
-    file = flist[ff]
+    #Get list of final images
+    namestr = '/Users/courtney/Documents/data/shaneAO/*/reduced*/*/*/contrast_curve.csv'
+    flist = glob.glob(namestr)
 
-    #Use current filename to set final filename
-    parts = file.split('/')
-    filt = parts[-2] #filter is second to last part of filename
-    tic = parts[-3] #TIC is third to last part of filename
-    night = parts[-4].split('reduced_')[1] #reduced_[NIGHT] is fourth to last part of filename
+    logger.info('Files: ', len(flist))
 
+    #Loop through these files!
+    counter = 0
+    for ff in np.arange(len(flist)):
 
-    #Don't include Kepler or K2 targets
-    if tic[0] == 'K': #Catch Kepler or K2 prefixes
-        continue
-    if tic[0] == 'E': #catch EPIC prefixes
-        continue
+        file = flist[ff]
 
-    #Remove T or TIC prefix
-    if 'TIC' in tic:
-        if verbose:
-            print('TIC name: ', tic)
-        tic = tic.split('TIC')[1]
-        if verbose:
-            print('renamed as : ', tic)
+        #Use current filename to set final filename
+        parts = file.split('/')
+        filt = parts[-2] #filter is second to last part of filename
+        tic = parts[-3] #TIC is third to last part of filename
+        night = parts[-4].split('reduced_')[1] #reduced_[NIGHT] is fourth to last part of filename
 
 
-    if 'T' in tic:
-        if verbose:
-            print('T name: ', tic)
-        tic = tic.split('T')[1]
-        if verbose:
-            print('renamed as : ', tic)
+        #Don't include Kepler or K2 targets
+        if tic[0] == 'K': #Catch Kepler or K2 prefixes
+            continue
+        if tic[0] == 'E': #catch EPIC prefixes
+            continue
 
-    #Recast to drop leading zeros and spaces
-    tic = str(int(tic))
+        #Remove T or TIC prefix
+        if 'TIC' in tic:
+            logging.info('TIC name: ', tic)
 
-    #Set output file
-    outname = tic+'_'+filt+'_'+night+'_contrast_curve.csv'
+            tic = tic.split('TIC')[1]
 
-    #Read in the contrast curve
-    c = pd.read_csv(file)
+            logging.info('renamed as : ', tic)
 
-    #Drop the error column
-    c = c[['arcsec','dmag']]
 
-    #Don't keep any rows with missing values
-    c = c.dropna()
 
-    #Write TRICERATOPS-friendly output file
-    c.to_csv(tridir+outname,index=False,header=False)
-    counter += 1
+        if 'T' in tic:
+            logging.info('T name: ', tic)
 
-print('Saved ', counter, ' contrast curves in TRICERATOPS format.')
+            tic = tic.split('T')[1]
+
+            logging.info('renamed as : ', tic)
+
+
+        #Recast to drop leading zeros and spaces
+        tic = str(int(tic))
+
+        #Set output file
+        outname = tic+'_'+filt+'_'+night+'_contrast_curve.csv'
+
+        #Read in the contrast curve
+        c = pd.read_csv(file)
+
+        #Drop the error column
+        c = c[['arcsec','dmag']]
+
+        #Don't keep any rows with missing values
+        c = c.dropna()
+
+        #Write TRICERATOPS-friendly output file
+        c.to_csv(tridir+outname,index=False,header=False)
+        counter += 1
+
+    logger.info('Saved '+str(counter)+' contrast curves in TRICERATOPS format.')
\ No newline at end of file
diff --git a/src/simmer/plotting.py b/src/simmer/plotting.py
index acc4c14..f0c4d2c 100644
--- a/src/simmer/plotting.py
+++ b/src/simmer/plotting.py
@@ -10,6 +10,9 @@
 
 from .schemas import read_yml as read
 
+import logging
+logger = logging.getLogger('simmer')
+
 plot_config = None
 
 
@@ -275,7 +278,7 @@ def plot_all_stars(func, snames=snames, filts=filts, nrows=4):
         fig, cim = plot_few(func, snames=snames)
 
     elif array_len > 50:
-        print("Too many images to plot.")
+        logger.info("Too many images to plot.")
         return
     else:  # 11 images? 13 images? make it 4xn
         fig, cim = plot_many(func, snames=snames)
diff --git a/src/simmer/registration.py b/src/simmer/registration.py
index 345ef7b..fd9734f 100644
--- a/src/simmer/registration.py
+++ b/src/simmer/registration.py
@@ -16,6 +16,9 @@
 
 from simmer.analyze_image import *
 
+import logging
+logger = logging.getLogger('simmer')
+
 
 def roll_shift(image, shifts, cval=0.0):
     """
@@ -365,9 +368,9 @@ def shift_bruteforce(image, base_position=None, max_shift=200, verbose=False):
     jlo = np.max([0, base_position[1] - max_shift])
     jhi  = np.min([imshape[0],base_position[1]+max_shift])
     if max_shift == 0:
-        if verbose:
-            print('ERROR: Max shiftset to 0. Considering full image.')
-            print('       Requested max_shift: ', max_shift)
+        logger.error('ERROR: Max shiftset to 0. Considering full image.')
+        logger.debug('       Requested max_shift: ', max_shift)
+
     else:
         masked_image= image.copy()*0.
         masked_image[ilo:ihi, jlo:jhi] = image[ilo:ihi, jlo:jhi]
diff --git a/src/simmer/run_night.py b/src/simmer/run_night.py
index fbced08..eabd8c1 100644
--- a/src/simmer/run_night.py
+++ b/src/simmer/run_night.py
@@ -11,8 +11,14 @@
 import simmer.check_logsheet as check
 import simmer.add_dark_exp as ad
 import simmer.create_config as config
+import logging
+
+logger = logging.getLogger('simmer')
 
 def run_night(wantdate, add_darks=True, just_images=False, sep_skies=False,  skip_reduction=False, verbose=False):
+    if verbose:
+        logger.setLevel(logging.DEBUG)
+
     #wantdate = desired night to reduce. Format is 'YYYY-MM-DD' (e.g., '2019-09-13')
 
     #Change these for your local installation
@@ -36,9 +42,8 @@ def run_night(wantdate, add_darks=True, just_images=False, sep_skies=False,  ski
 
     #Check that logsheet exists
     if os.path.isfile(log_name) == False:
-        print('ERROR. Logsheet does not exist. Attempted to find ', log_name, '.')
-    if verbose == True:
-        print('Reading logsheet ', log_name)
+        logger.error('ERROR. Logsheet does not exist. Attempted to find ', log_name, '.')
+    logger.debug('Reading logsheet ', log_name)
 
     #Set up config file from logsheet
     tab = None # this is a CSV, so we don't have a tab name of interest.
@@ -47,8 +52,8 @@ def run_night(wantdate, add_darks=True, just_images=False, sep_skies=False,  ski
     #add darks
     if add_darks == True:
         log_name_with_darks = ad.add_dark_exp(inst, log_name, rawdir, tab,)
-        if verbose == True:
-            print('Saving logsheet with darks added as: ', log_name_with_darks)
+        logger.debug('Saving logsheet with darks added as: ', log_name_with_darks)
+
     else:
         log_name_with_darks = log_name
 
@@ -58,7 +63,7 @@ def run_night(wantdate, add_darks=True, just_images=False, sep_skies=False,  ski
 
     #Reduce the data!
     if skip_reduction == True:
-        print('files exist')
+        logger.info('Skipping reduction')
         return config_file
     else:
         drivers.all_driver(inst, config_file, rawdir, reddir, just_images=just_images, sep_skies=sep_skies, verbose=verbose)
diff --git a/src/simmer/search_headers.py b/src/simmer/search_headers.py
index a809649..c170516 100644
--- a/src/simmer/search_headers.py
+++ b/src/simmer/search_headers.py
@@ -7,6 +7,8 @@
 from glob import glob
 
 import astropy.io.fits as pyfits
+import logging
+logger = logging.getLogger('simmer')
 
 
 def search_headers(raw_dir, write_dir=None):
@@ -50,7 +52,7 @@ def search_headers(raw_dir, write_dir=None):
                 textfile.write(f"{file}\n")
 
         else:
-            print(f"Header Incomplete in {file}!!! ")
+            logger.error(f"Header Incomplete in {file}!!! ")
             textfile.write(f"{file} HEADER INCOMPLETE\n")
 
     textfile.close()