######All saved in SI units###############################
import numpy as np
#constants
c = 299792458.
kB = 1.3806504e-23
h = 6.62606957e-34
hbar = h /(2.0*np.pi)
G = 6.67428e-11
#angles
radian = 1.0
degree = np.pi/180.
arcmin = degree/60.
arcsec = arcmin/60.
#mass
PlanckMass = np.sqrt(hbar*c/G)
ReducedPlanckMass = PlanckMass/np.sqrt(8*np.pi)
#length
PlanckLength = np.sqrt(G*hbar/c**3)
pc = 3.08568025e16
kpc = 1.e3*pc
Mpc = 1.e6*pc
Gpc = 1.e9*pc
#time
PlanckTime =  PlanckLength/c #np.sqrt(G*hbar/c)/c**2
day = 3600.*24.
year = day*365.2425
kyr = 1.e3*year
Myr = 1.e6*year
Gyr = 1.e9*year
#energy
PlanckEnergy = PlanckMass*c**2 #np.sqrt(hbar*c/G)*c**2
eV = 1.60217648740e-19
keV = 1.e3 * eV
MeV = 1.e6 * eV
GeV = 1.e9 * eV
#temprature
PlanckTemperature = PlanckEnergy/kB #np.sqrt(hbar*c/G)*c**2/k_B
#Atomic constants
AtomicMassUnit = 1.66053878283e-27
ElectronMass = AtomicMassUnit/1822.8884845
ProtonMass = 1.00727646681290 * AtomicMassUnit
NeutronMass = 1.0086649160043 * AtomicMassUnit
HydrogenBindingEnergy = 13.605698 * eV
HydrogenMass = ElectronMass + ProtonMass - HydrogenBindingEnergy/c**2
HeliumMass = 3.9715*HydrogenMass
#Astronomical constants
LightYear = c*year
SunMass = 1.98892e30
EarthMass = 5.97219e24
H0Unit = 1.e3/Mpc  #km/s/Mpc
#################################################
print("PlanckTime = ", PlanckTime)
print("PlanckLength = ", PlanckLength)
print("PlanckMass = ", PlanckMass)
print("PlanckEnergy = ", PlanckEnergy)
print("PlanckTemperature = ", PlanckTemperature)
##=========User block, write your code here ========
print(r'Mpc^{-1} = ', (PlanckLength/Mpc)/(eV/PlanckEnergy), ' eV')

H0 = 70 * H0Unit * PlanckTime ##in Planck units
TCMB = 2.73/PlanckTemperature  ## in Planck units
rho_gamma = (np.pi**2/15)*TCMB**4 ##in Planck units
print(r'Omega_gamma = ', (8 * np.pi * rho_gamma)/(3.0*H0**2))
#===============End ================================