meanTempChangePerOrbit = np.mean(temps)

morbits = 1e9
orbitJump = 5e2

orbits = np.arange(0, morbits, orbitJump)

hotSideTemp = np.zeros_like(orbits)
coldSideTemp = np.zeros_like(orbits)

hottemp = BDTEMP
coldtemp = BDTEMP

period = 6*3600 #seconds
longTimes = np.linspace(0, morbits*period, orbits.shape[0])

bdmass = 0.075 * SMASS
intrinsicCoolingCoef = 0.9
initialCooling = hottemp - cool_blackbody(BDRAD, hottemp, bdmass, period*orbitJump, 0)
internalEnergy = (3/2) * N * kb * initialCooling * intrinsicCoolingCoef


for oID, orbit in tqdm(enumerate(orbits), total=orbits.shape[0]):
    hottemp += (orbitID + 1) * meanTempChangePerOrbit
    hottemp = cool_blackbody(BDRAD, hottemp, bdmass, orbitJump * period, internalEnergy) 
    hotSideTemp[oID] = hottemp