def calibrate_2band(instr1, instr2, airmass1, airmass2, coeff1, coeff2,
        zero_key='zero', color_key='color', extinct_key='extinction'):
    """
    This solves the set of equations:
        i_0 = i + A_i + C_i(i-z) + k_i X
        z_0 = z + A_z + C_z(z-i) + k_z X
    where i_0 and z_0 are the instrumental magnitudes, A_i and A_z are the 
    zero points, C_i and C_z are the color terms, k_i and k_z are the 
    atmospheric coefficients, and X is the airmass.
    
    The solution is of the form:
        (1+C_i)i = b_i + C_i z
        (1+C_z)z = b_z + C_z i
    where
        b_i = i_0 - A_i - k_i X
        b_z = z_0 - A_z - k_z X
    so that
        i = (C_i b_z + C_z b_i + b_i) / d
        z = (C_z b_i + C_i b_z + b_z) / d
    where
        d = (1+C_i+C_z)
    
    Parameters
    ----------
    instr1: array-like
        Instrumental magnitudes of filter 1
    instr2: array-like
        Instrumental magnitudes of filter 2
    airmass1: array-like
        Airmass for each observation in filter 1
    airmass2: array-like
        Airmass for each observation in filter 2
    coeff1: array-like
        List of coeffients for calibrating instrumental magnitudes for 
        instrument 1.
            * coeff1[0]: zeropoint
            * coeff1[1]: color coeffcient
            * coeff1[2]: extinction coefficient
    coeff2: array-like
        List of coeffients for calibrating instrumental magnitudes for 
        instrument 2
    
    returns
    -------
    mag1: array-like
        Calibrated magnitude 1
    mag2: array-like
        Calibrated magnitude 2
    """
    b1 = instr1 - coeff1[zero_key] - coeff1[extinct_key]*airmass1
    b2 = instr2 - coeff2[zero_key] - coeff2[extinct_key]*airmass2
    d = 1 + coeff1[color_key] + coeff2[color_key]
    mag1 = (coeff1[color_key]*b2 + b1*(1+coeff2[color_key])) / d
    mag2 = (coeff2[color_key]*b1 + b2*(1+coeff1[color_key])) / d
    return (mag1,mag2)