#!/usr/bin/env python

import os
import sys
import pyfits
import numpy


if __name__ == "__main__":

    blue_list = sys.argv[1]
    red_list = sys.argv[2]

    color_list = sys.argv[3]


    with open(blue_list, "r") as f:
        blue_fn_list = [l.split()[0] for l in f.readlines()]
    with open(red_list, "r") as f:
        red_fn_list = [l.split()[0] for l in f.readlines()]
    with open(color_list, "r") as f:
        color_fn_list = [l.split()[0] for l in f.readlines()]

    for idx in range(len(blue_fn_list)):
        print "Working on file %d of %d" % (idx+1, len(blue_fn_list))

        blue_fn = blue_fn_list[idx]
        red_fn = red_fn_list[idx]
        color_fn = color_fn_list[idx]

        print "opening", blue_fn
        blue_hdu = pyfits.open(blue_fn)
        print "opening", red_fn
        red_hdu = pyfits.open(red_fn)

        zp = 27
        print "computing blue magnitudes"
        blue = -2.5*numpy.log10(blue_hdu['IMAGE'].data) + zp
        print "computing red magnitudes"
        red = -2.5*numpy.log10(red_hdu['IMAGE'].data) + zp


        bad = (blue_hdu['IMAGE'].data <= 0) | (red_hdu['IMAGE'].data <= 0)

        print "computing color"
        # color = red - blue
        color = blue - red
        color[bad] = numpy.NaN

        print "writing output to", color_fn
        pyfits.PrimaryHDU(header=blue_hdu['IMAGE'].header, data=color).writeto(color_fn, clobber=True)

        print "\n\n"

    print "all done!"