def METHOD_NAME(day_imgfile, night_imgfile, sublon, sublat, out_imgfile):
    """Given filenames for a day and a night image,
       and a subsolar longitude and latitude,
       combine the two images and write the result to a new image file.
       But if the image exists already, skip it, don't overwrite.
    """
    if os.path.exists(out_imgfile):
        return
    dayimg = Image.open(day_imgfile)
    nightimg = Image.open(night_imgfile)
    # *img.shape is height, width, depth
    w, h = dayimg.size
    halfwidth = w/2
    halfheight = h/2
    lon_conversion = 180. / halfwidth
    lat_conversion = 90. / halfheight
    # Iterate over pixels
    dayimg_px = dayimg.load()
    nightimg_px = nightimg.load()
    for y in range(h):
        lat = (halfheight - y) * lat_conversion
        for x in range(w):
            # pxl is the RGB 3-element array at (x, y)
            lon = (x - halfwidth) * lon_conversion
            # The Haversine angle will be between 0 and 180 degrees.
            havangle = haversine_angle(sublon, sublat, lon, lat)
            if havangle > 90 + BLURANGLE:
                # full night
                dayimg_px[x, y] = nightimg_px[x, y]
            elif havangle > 90 - BLURANGLE and havangle < 90 + BLURANGLE:
                # twilight, blur zone.
                blur = (90 + BLURANGLE - havangle) / BLURANGLE / 2.
                # Each pixel is a tuple, so make a new tuple:
                newpx = []
                for i, c in enumerate(dayimg_px[x, y]):
                    newpx.append(int(dayimg_px[x, y][i] * blur
                                     + nightimg_px[x, y][i] * (1. - blur)))
                dayimg_px[x, y] = tuple(newpx)
            # else it's daylight, don't do annything
    # dayimg.show()
    dayimg.save(out_imgfile)