import numpy
from particles import Particle
from evolution import Universe
from library.SM import particles as SMP
from common import Params, UNITS
T_final = 100 * UNITS.keV
params = Params(T=100 * UNITS.MeV,
dx=1e-2 * UNITS.MeV)
universe = Universe(params=params, logfile="tests/photon_universe/log.txt")
photon = Particle(**SMP.photon)
universe.add_particles(photon)
universe.evolve(T_final)
initial_aT = universe.data['aT'][0]
print "a * T is conserved: {}".format(all([initial_aT == value for value in universe.data['aT']]))
initial_a = universe.data['a'][len(universe.data['a'])/2]
initial_t = universe.data['t'][len(universe.data['a'])/2] / UNITS.s
last_a = universe.data['a'][-1]
last_t = universe.data['t'][-1] / UNITS.s
print "a scaling discrepancy is: {:.2f}%"\
.format(100 * (last_a / initial_a / numpy.sqrt(last_t / initial_t) - 1))
universe.graphics.save(__file__)
raw_input("...")