def cool_blackbody(radius, initialTemperature, mass, delta_time, internal_energy):
    N = estimateParticlesBBN(mass*BDATMMFRAC) # assume all energy is distributed through the atmosphere
    Ui = 3/2 * N * kb * initialTemperature
    
    initial_luminosity = 4 * np.pi * radius**2 * sb * (initialTemperature)**4
    energy_lost = initial_luminosity * delta_time

    Uf = (Ui - energy_lost) + internal_energy
    finalTemperature = (2/3 * Uf)/(N*kb)
    return finalTemperature