function bplanck(wav, temp)
        implicit double precision (o-z, a-h)
        parameter (one = 1.d0, c1 = 1.192d-5, c2 = 1.439565d0)
        c3 = dexp(- ((c2 * wav) / temp))
        bplanck = ((c1 * (wav ** 3)) * c3) / (one - c3)
        return 
        end