#!/usr/bin/env python

import os
import sys
import pyfits
import numpy


if __name__ == "__main__":

    blue_fn = sys.argv[1]
    red_fn = sys.argv[2]

    color_fn = sys.argv[3]

    blue_hdu = pyfits.open(blue_fn)
    red_hdu = pyfits.open(red_fn)

    zp = 27
    print "computing blue magnitudes"
    blue = -2.5*numpy.log10(blue_hdu[0].data) + zp
    print "computing red magnitudes"
    red = -2.5*numpy.log10(red_hdu[0].data) + zp
    bad = (blue_hdu[0].data <= 0) | (red_hdu[0].data <= 0)
    
    print "computing color"
    # color = red - blue
    color = blue - red
    color[bad] = numpy.NaN

    print "writing output to", color_fn
    pyfits.PrimaryHDU(data=color).writeto(color_fn, clobber=True)


    print "all done!"