30:T5f8,md`## Formalism Background quantities are derived from: ${tex.block` E(z) = \left[ \Omega_{\text{DE},0} f_\Lambda(z) + \Omega_{k,0} (1+z)^2 + \Omega_{m,0} (1+z)^3 + \Omega_{\gamma,0} (1 + f_\nu(z)) (1+z)^4 \right]^{1/2} `} Specifically, line-of-sight comoving distance is: ${tex.block` d(z_1,z_2) = \frac{c}{H_0} I_0(z_1,z_2) `} and the corresponding lookback time is: ${tex.block` t(z_1,z_2) = \frac{1}{H_0} I_1(z_1,z_2) `} with ${tex.block` I_p(z) \equiv \int_{z_1}^{z_2} \frac{dz}{(1+z)^p E(z)} = \int_{1/(1+z_2)}^{1/(1+z_1)} \frac{da}{a^{2-p} E(a)} \; . `} The particle horizon distance and age of the universe are the limits of these quantitites as ${tex`z_2\rightarrow\infty`} or ${tex`a_1\rightarrow 0`}. Transverse comoving distances are related to the line-of-sight comoving distance above via: ${tex.block` S_\kappa(r) = \begin{cases} R_0 \sin\left(r/R_0 \right) & \kappa > 0 \\ r & \kappa = 0 \\ R_0 \sinh\left(r/R_0 \right) & \kappa < 0 \\ \end{cases} `} with ${tex.block` R_0 \equiv \frac{c}{H_0} \left|\Omega_{k,0}\right|^{-1/2} \quad, \quad \kappa = \text{sign}\left(\Omega_{\text{tot},0} - 1\right) = \text{sign}\left(-\Omega_{k,0}\right) \; . `} Dark energy density scales according to: ${tex.block` \Omega_\text{DE}(z) = \Omega_{\text{DE},0} f_\Lambda(z) `} with ${tex.block` f_\text{DE}(z) \equiv \exp\left( 3 \int_a^1 \frac{da'}{a'} \left[ 1 + w(a')\right] \right) =\begin{cases} 1 & w(z)=-1 \\ (1+z)^{3(1+w_0+w_a)} \exp\left(\frac{-3 w_a z}{1+z}\right) & w(z) = w_0 + \frac{w_a z}{1+z} \end{cases} `} `31:T526,class LambdaCDM extends FLRW { // Concrete subclass with constant w(z)=-1 and arbitrary curvature. Must specify Ok0. // If w0 or wa are present in the ctor params, they must have the values w0=-1, wa=0. constructor(params={}) { if((params.w0 !== undefined) && (params.w0 != -1)) throw Error('Cannot specify w0 != -1 for LambdaCDM'); if((params.wa !== undefined) && (params.wa != 0)) throw Error('Cannot specify wa != 0 for LambdaCDM'); super(params); } w(z) { return -1; } de_density_scale(z) { return 1; } de_density_scale_a4(a) { return a**4; } // Return a model with the same early-time evolution a(t) as the current model, // but no dark energy, which is useful for demonstrating the accelerated // expansion due to dark energy. non_accelerated() { const r = 1 - this.Ode0; // Scale H0 and the Omegas to maintain the same physical energy densities associated // with matter and curvature. No change to Tcmb since Tcmb**4 sets the physical energy // density, which we want to maintain. const params = { H0: this._H0 * Math.sqrt(r), Om0: this.Om0 / r, Ob0: this.Ob0 / r, Ok0: this.Ok0 / r, Tcmb0: this._Tcmb0, Neff: this.Neff, m_nu: this.m_nu, }; return new LambdaCDM(params); } }32:T1628,class FLRW extends BaseModel { // This is an abstract base class with an unspecified dark energy sector. // Concrete subclasses must implement w(z), de_energy_scale(z), de_energy_scale_a4(a). constructor({H0,Om0,Ob0,Ok0,Tcmb0,Neff,m_nu=0}={}) { super({H0:H0,Ok0:Ok0}) if(new.target === FLRW) throw Error('Cannot instantiate abstract base class FLRW.'); if(this.w === undefined) throw Error('Subclass does not implement w(z)'); if(this.de_density_scale === undefined) throw Error('Subclass does not implement de_density_scale(z)'); if(this.de_density_scale_a4 === undefined) throw Error('Subclass does not implement de_density_scale_a4(a)'); // == Initialize the matter sector. this._Om0 = Om0; if(Ob0 > Om0) throw Error('Expected Ob0 <= Om0.'); this._Ob0 = Ob0; // == Initialize the photon sector. this._Tcmb0 = Tcmb0; // Kelvin // Compute the present photon relative density. this._Ogamma0 = Const.planck_alpha * this._Tcmb0**4 / this._critical_density0; // == Initialize the neutrino sector. // Compute the corresponding present-day neutrino temperature. this._Tnu0 = (4/11)**(1/3) * this._Tcmb0; // Kelvin if(Neff < 0) throw Error('Expected Neff >= 0.'); this._Neff = Neff; this._N_nu = Math.floor(Neff); if(typeof m_nu === 'number') { this._m_nu = new Array(this._N_nu).fill(m_nu); // eV } else if(!Array.isArray(m_nu)) { throw Error('Expected scalar or array for m_nu.'); } else if(m_nu.length != this._N_nu) { throw Error(`Expected ${this._N_nu} values for m_nu.`); } else { // should check that each m_nu is a number and >= 0 here... this._m_nu = m_nu; // eV } if(this._N_nu > 0) { const kB = 8.617333262145179e-05; // in eV K this._nu_y = this._m_nu.map(m => m / (kB * this._Tnu0)); // divide by 1+z to get arg to f(y) // Fitting formula from eqn (26) of Komatsu 2011, i.e. f(y) = (1 + (A y)**p) ** (1/p) // Note that we use A = 0.3173 to match astropy.cosmology.FLRW instead of the more // precise value A = 180 zeta(3) / (7 pi**4) in the paper. this._nu_f = y => (1 + (0.3173 * y) ** 1.83) ** 0.54644808743; const nu_const = 7/8 * (4/11)**(4/3) * this._Neff / this._N_nu; this.nu_relative_density = z => { let answer = 0; this._nu_y.forEach(y => { answer += this._nu_f(y / (1 + z)); }); return nu_const * answer; }; } else { this.nu_relative_density = z => 0; } this._Onu0 = this.nu_relative_density(0) * this._Ogamma0; // == Initialize the dark-energy sector. // Compute the present dark energy relative density. this._Ode0 = 1 - this.Ok0 - this._Om0 - this._Ogamma0 - this._Onu0; } // Getters for present-day (constant) dimensionless parameter values. get Om0() { return this._Om0; } get Ob0() { return this._Ob0; } get Ode0() { return this._Ode0; } get Ogamma0() { return this._Ogamma0; } get Onu0() { return this._Onu0; } get Neff() { return this._Neff; } // Memoized getters for constant quantities with dimensions. get Tcmb0() { return this._Tcmb0q = this._Tcmb0q || quantity(this._Tcmb0, 'K'); } get Tnu0() { return this._Tnu0q = this._Tnu0q || quantity(this._Tnu0, 'K'); } get critical_density0() { return this._critical_density0q = this._critical_density0q || quantity(this._critical_density0, 'g/cm^3'); } // Methods to compute parameters that evolve as some power of E(z). Om(z) { return this._Om0 * (1+z)**3 / this.efuncsq(z); } Odm(z) { return (this._Om0 - this._Ob0) * (1+z)**3 / this.efuncsq(z); } Ob(z) { return this._Ob0 * (1+z)**3 / this.efuncsq(z); } Ode(z) { return this._Ode0 * this.de_density_scale(z) / this.efuncsq(z); } Ogamma(z) { return this._Ogamma0 * (1+z)**4 / this.efuncsq(z); } Onu(z) { return this.Ogamma(z) * this.nu_relative_density(z); } Tgamma(z) { return this._Tcmb0 * (1 + z); } Tnu(z) { return this._Tnu0 * (1 + z); } ngamma(z) { return Const.planck_beta * this.Tgamma(z) ** 3; } // photons/cm^3 // Convert baryon and DM energy densities to equivalent number densities assuming non-relativistic // particles of mass mrel relative to the proton mass. nb(z, mrel=1) { return this.Ob(z) * this.critical_density(z) / (mrel * Const.proton_mass) }; // baryons/cm^3 ndm(z, mrel=1) { return this.Odm(z) * this.critical_density(z) / (mrel * Const.proton_mass) }; // DM/cm^3 efuncsq(z) { // Define the evolution factor E(z)**2. const Or = this._Ogamma0 * (1 + this.nu_relative_density(z)); const zp1 = 1 + z; return this._Ode0 * this.de_density_scale(z) + zp1 ** 2 * (this._Ok0 + zp1 * (this._Om0 + zp1 * Or)); } integrand_z(z, p) { // Evaluate 1/[(1+z)^p E(z)]. Use p=1 to evaluate lookback time in units of the Hubble time // or p=0 to evaluate comoving distance in units of the Hubble distance. Use q=2-p for // the corresponding integral of integrand_a(). const Or = this._Ogamma0 * (1 + this.nu_relative_density(z)); const zp1 = 1 + z; return 1 / (zp1**p * Math.sqrt( this._Ode0 * this.de_density_scale(z) + zp1 ** 2 * (this._Ok0 + zp1 * (this._Om0 + zp1 * Or)))); } integrand_a(a, q) { // Evaluate 1/[a^q E(a)]. Use q=1 to evaluate lookback time in units of the Hubble time // or q=2 to evalute comoving distance in units of the Hubble distance. Use p=2-q for // the corresponding integral of integrand_z(). const Or = this._Ogamma0 * (1 + this.nu_relative_density(1/a - 1)); return 1 / (a**(q-2) * Math.sqrt( this._Ode0 * this.de_density_scale_a4(a) + Or + a * (this._Om0 + a * this._Ok0))); } }33:T151e,class FluidsModel extends BaseModel { // Implement a model with arbitrary radiation, matter and dark energy // fluid components with fixed 1/a**(3*(1+w)) evolution. Compared with // FLRW, additional evolution due to non-constant w for dark energy // or massive neutrinos is not included, but models with no "big bang" // (a=0) or where adot=0, leading to big crunch/bounce evolution, can // be explored via the solve method. Note that the BaseModel age // will only be finite if at least one of (Or0,Om0,Ok0) is non zero. // Similarly, a finite horizon size requires Or0 nonzero. constructor(Or0, Om0, Ode0, H0) { const Ok0 = 1 - (Or0 + Om0 + Ode0); super({H0:H0, Ok0:Ok0}) this.Or0 = Or0; this.Om0 = Om0; this.Ode0 = Ode0; //this.Omega = new Float64Array([Or0, Om0, Ok0, Ode0]); //this.pow = new Float64Array([-4, -3, -2, 0]); // Find any positive scales where adot = 0. this.aminmax = realsolns(quarticsolve(Ode0, 0, Ok0, Om0, Or0)).filter(a => a > 1e-10); const conversion = Const.sec_to_Gyr / (H0 * Const.H0units_to_invs); this.convert = u => conversion * u; //this.Esq = a => a==0 ? 0 : Math.max(1e-8, Ode0 + (Ok0 + (Om0 + Or0 / a) / a) / (a * a)); this._integrand = function(a) { if(a == 0) return 0; const a2 = a * a; return 1 / Math.sqrt(Or0 / a2 + Om0 / a + Ok0 + Ode0 * a2); } // Evaluate the slope of the polynomial whose roots correspond to points where adot=0. this.polyslope = a => Om0 + a * (2 * Ok0 + 4 * Ode0 * a * a); } efuncsq(z) { const zp1 = 1 + z; return this.Ode0 + zp1*zp1 * (this.Ok0 + zp1 * (this.Om0 + zp1 * this.Or0)); } integrand_z(z, p) { // Evaluate 1/[(1+z)^p E(z)]. Use p=1 to evaluate lookback time in units of the Hubble time // or p=0 to evaluate comoving distance in units of the Hubble distance. Use q=2-p for // the corresponding integral of integrand_a(). const zp1 = 1 + z; return 1 / (zp1**p * Math.sqrt(this.Ode0 + zp1 ** 2 * (this.Ok0 + zp1 * (this.Om0 + zp1 * this.Or0)))); } integrand_a(a, q) { // Evaluate 1/[a^q E(a)]. Use q=1 to evaluate lookback time in units of the Hubble time // or q=2 to evalute comoving distance in units of the Hubble distance. Use p=2-q for // the corresponding integral of integrand_z(). return this.Or0 == 0 ? (this.Om0 == 0 ? 1 / (a**(q-1) * Math.sqrt(this.Ok0 + a*a * this.Ode0)) : 1 / (a**(q-1.5) * Math.sqrt(this.Om0 + a * (this.Ok0 + a*a * this.Ode0)))) : 1 / (a**(q-2) * Math.sqrt(this.Or0 + a * (this.Om0 + a * (this.Ok0 + a*a * this.Ode0)))); } _integrate(agrid, nquad) { return Float64Array.from( {length:agrid.length}, (_,i) => i==0 ? 0 : quad(this._integrand, agrid[i-1], agrid[i], nquad)); } divergent(a0, da) { // Use series approx for analytic integral of da/a over [a0-da,a0] or [a0,a0+da] // where the integrand diverges at a0. const slope0 = this.polyslope(a0); const sign = Math.sign(slope0); // +1 is [a0,a0+da], -1 is [a0-da,a0]. return 2/3 * da * (3 * a0 + sign * da) / Math.sqrt(sign * da * slope0); } solve({amin=0, amax=2, ngrid=32, nquad=16, reflect=true}={}) { if(amin > 1 || amax < 1) throw new Error('[amin,amax] must bracket 1.'); // Trim the range if there are min/max points. const a1cut = Math.max(amin-1, Math.max(...this.aminmax.filter(a => a < 1))); const a2cut = Math.min(amax+1, Math.min(...this.aminmax.filter(a => a > 1))); // Integrate over [a1,1]. const a1grid = linspace(Math.max(amin,a1cut), 1, ngrid, {dtype:Array}); let t1grid = this._integrate(a1grid, nquad); if(a1cut >= amin) { // Solution hits adot=0 so replace the first interval. t1grid[1] = this.divergent(a1cut, a1grid[1] - a1grid[0]); } // Adjust times to be relative to a=1. t1grid = cumsum(t1grid); const t0 = t1grid[ngrid - 1]; t1grid = Array.from(t1grid).slice(0, -1).map(t => t - t0); // Integrate over [1,a2]. const a2grid = linspace(1, Math.min(amax,a2cut), ngrid); let t2grid = this._integrate(a2grid, nquad); if(a2cut <= amax) { // Solution hits adot=0 so replace the last interval. t2grid[ngrid-1] = this.divergent(a2cut, a2grid[ngrid-1] - a2grid[ngrid-2]); } t2grid = cumsum(t2grid); // Assemble the two halves. this.agrid = arraycat(a1grid.slice(0, -1), a2grid); this.tgrid = arraycat(t1grid, t2grid).map(this.convert); this.agrid_unreflected = Array.from(this.agrid); this.tgrid_unreflected = Array.from(this.tgrid); if(reflect & ((a1cut >= amin) || (a2cut <= amax))) { if((a1cut >= amin) && (a2cut <= amax)) { console.log('Can this ever happen?'); } else if(a1cut >= amin) { // Reflect the solution from [a1cut, amax] about a1cut. this.agrid = arraycat(Array.from(this.agrid).slice(1).reverse(), this.agrid); const t0 = this.tgrid[0]; this.tgrid = arraycat(Array.from(this.tgrid).slice(1).reverse().map(t => 2*t0-t), this.tgrid); } else if(a2cut <= amax) { // Reflect the solution from [amin, a2cut] about a2cut. this.agrid = arraycat(this.agrid, Array.from(this.agrid).slice(1).reverse()); const t0 = this.tgrid[2*ngrid-2]; this.tgrid = arraycat(this.tgrid, Array.from(this.tgrid).slice(1).reverse().map(t => 2*t0-t)); } } return this; } }34:T18e5,class BaseModel { // This an abstract base class for a model with unspecified components. // Concrete subclasses must implement efuncsq, integrand_a, integrand_z. constructor({H0,Ok0}={}) { this._H0 = H0; // (km/s)/Mpc // ...present Hubble distance c/H0 in Mpc and Hubble time 1/H0 in Gyr. this._hubble_distance = Const.c_kms / this._H0; this._hubble_time = Const.sec_to_Gyr / (Const.H0units_to_invs * this._H0); // ...present critical density in g/cm^3. this._critical_density0 = Const.critdens_const * (this._H0 * Const.H0units_to_invs) ** 2; this._Ok0 = Ok0; // Curvature function for transverse comoving distances. this.kappa = Math.sign(-Ok0); // Omega_{tot,0} - 1 = -Ok0 this.R0 = this._hubble_distance / Math.sqrt(Math.abs(Ok0)); if(this.kappa > 0) { this.Sk = r => this.R0 * Math.sin(r / this.R0); } else if(this.kappa < 0) { this.Sk = r => this.R0 * Math.sinh(r / this.R0); } else { // this.kappa == 0 this.Sk = r => r; } // Use introspection to identify derived quantities that involve an integral over z // so are written as functions of (z1,z2) or else are direct functions of (z). // These regexps could be modified to allow additional optional args, but then // be careful since "(z1,z2," starts with "(z". const is_integrated_regexp = /$z1,\s*z2$/; this.is_integrated = f => is_integrated_regexp.test(f.toString()); const is_direct_regexp = /$z$/; this.is_direct = f => is_direct_regexp.test(f.toString()); } // Getters for present-day (constant) dimensionless parameter values. get h() { return this._H0 / 100; } get Ok0() { return this._Ok0; } // Memoized getters for constant quantities with dimensions. get H0() { return this._H0q = this._H0q || quantity(this._H0, '(km/s)/Mpc'); } get hubble_distance() { return this._hubble_distanceq = this._hubble_distanceq || quantity(this._hubble_distance, 'Mpc'); } get hubble_time() { return this._hubble_timeq = this._hubble_timeq || quantity(this._hubble_time, 'Gyr'); } // Methods to compute parameters that evolve as some power of E(z). H(z) { return this._H0 * Math.sqrt(this.efuncsq(z)); } // Calculate adot(z) in inverse Gyr adot(z) { return this.H(z) / (1 + z) * Const.H0units_to_invs / Const.sec_to_Gyr; } Ok(z) { return this._Ok0 * (1+z)**2 / this.efuncsq(z); } critical_density(z) { // Result is in g/cm^3. return this._critical_density0 * this.efuncsq(z); } integral(p, z1, z2, zbreak=1) { // Estimate the integral of dz/[(1+z)^p E(z)] over [z1,z2]. // Uses the equivalent integral of da/[a^(2-p) E(a)] over [1/(1+z2),1/(1+z1)] // for any part of the integral above zbreak, for better convergence. // Uses Clenshaw-Curtis integration from @dkirkby/integrate. // The integration order is based on some simple heuristics that could be improved. // Swapping z1,z2 flips the sign of the result. const sign = Math.sign(z2-z1); if(sign == 0) return 0; if(sign < 0) [z1,z2] = [z2,z1]; // require z1 1) ? 64 : ((dz > 0.1) ? 32 : ((dz > 0.01) ? 16 : 8))); return sign * quad(z => this.integrand_z(z, p), z1, z2, nquad); } else if(z1 >= zbreak) { // Integrate over a. const a1 = 1/(1+z2), a2 = 1/(1+z1); const da = a2 - a1; const nquad = (a1 < 0.02) ? 128 : ((da > 0.5) ? 64 : ((da > 0.1) ? 32 : ((da > 0.01) ? 16 : 8))); return sign * quad(a => this.integrand_a(a, 2-p), a1, a2, nquad); } else { // Split into z and a integrals at zbreak. return this.integral(p, z1, zbreak) + this.integral(p, zbreak, z2); } } get age() { // Compute lookback time from a=1 to a=0 in Gyr. return this._age = this._age || this._hubble_time * this.integral(1, 0, Infinity); } get horizon() { // Compute the particle horizon distance in Mpc. return this._horizon = this._horizon || this._hubble_distance * this.integral(0, 0, Infinity); } comoving_distance(z1, z2) { // Calculate comoving distance in Mpc from z1 to z2. // Distances from z=0 to z>0 are >0 by convention. // If z1,z2 are swapped the sign of the result is flipped. return this._hubble_distance * this.integral(0, z1, z2); } comoving_transverse_distance(z1, z2) { // Calculate comoving transverse distance in Mpc from z1 to z2. return this.Sk(this.comoving_distance(z1, z2)); } angular_diameter_distance(z) { return this.comoving_transverse_distance(0, z) / (1 + z); } luminosity_distance(z) { return this.comoving_transverse_distance(0, z) * (1 + z); } lookback_time(z1, z2) { // Calculate lookback time in Gyr from z1 to z2. // Times from z=0 to z>0 are >0 by convention. // If z1,z2 are swapped the sign of the result is flipped. return this._hubble_time * this.integral(1, z1, z2); } lookback_distance(z1, z2) { // Lookback distance is c * lookback time in Gpc. return Const.c_Mpc_per_Gyr * this.lookback_time(z1, z2); } gridz(z, fname) { // Return an array of fname values corresponding to the input array of redshifts z. // Handles both integrated quantities (comoving_distance, etc) and quantities that // scale as some power of E(z). if(typeof this[fname] !== 'function') throw Error('Invalid method name: "' + fname + '".'); const method = this[fname]; const f = method.bind(this); if(this.is_direct(method)) { return z.map(z => f(z)); } else if(this.is_integrated(method)) { // Verify that z is strictly increasing. if(!z.every((_,i,z) => i==0 || z[i] > z[i-1])) throw Error('Values of z must be increasing.'); // Estimate the starting point. const f1 = f(0, z[0]); // Estimate separate integrals over each step. const df = Float64Array.from({length: z.length}, (_,i) => i==0? f1 : f(z[i-1], z[i])); // Return the cummulative sum relative to z[0]. return cumsum(df); } else throw Error('Method not a function of (z) or (z1,z2): "' + fname + '".'); } tidyz(z, fnames) { const data = {z:z}; for(const fname of fnames.split(/,\s*/)) data[fname] = this.gridz(z, fname); return tidy(data); } }35:T89b,Const = { class MemoizedConstants { // Memoized getters for useful unit conversion constants. get H0units_to_invs() { // Conversion from (km/s)/Mpc to 1/s. return this._H0units_to_invs = this._H0units_to_invs || quantity('(km/s)/Mpc', 's^-1').value; } get c_kms() { // Speed of light in km/s. return this._c_kms = this._c_kms || quantity('speedOfLight', 'km/s').value; } get c_Mpc_per_Gyr() { return this._c_Mpc_per_Gyr = this._c_Mpc_per_Gyr || quantity('speedOfLight', 'Mpc/Gyr').value; } get sec_to_Gyr() { // Conversion from s to Gyr. return this._sec_to_Gyr = this._sec_to_Gyr || quantity('1 s', 'Gyr').value; } get critdens_const() { // 3 / (8 pi G) in g s^2 / cm^3, which is useful for calculating the critical density. return this._critdens_const = this._critdens_const || quantity('3/(8 pi gravitationConstant)', 'g s^2/cm^3').value; } get planck_alpha() { // Return the constant of proportionality between the effective density (g/cm^3) of a black-body // radiation field (i.e. energy density / c^2) // at temperature T in Kelvin and the value of T^4, equal to pi^2/15 kB^4/(hbar^3 c^5). return this._planck_alpha = this._planck_alpha || quantity('(pi^2/15) boltzmann^4/(reducedPlanckConstant^3 speedOfLight^5)', 'g/(cm^3 K^4)').value; } get planck_beta() { // Return the constant of proportionality between the number density (/cm^3) of a black-body // radiation field at temperature T in Kelvin and the value of T^3, equal to // 2 zeta3/pi^2 kB^3/(hbar^3 c^3) with zeta3 = 1.2020569031595942854. return this._planck_beta = this._planck_beta || 1.2020569031595942854 * quantity( '(2/pi^2) boltzmann^3/(reducedPlanckConstant^3 speedOfLight^3)', 'cm^-3 K^-3').value; } get proton_mass() { // Proton mass in g. return this._proton_mass = this._proton_mass || quantity('protonMass', 'g').value; } get kB_evK() { // The Boltzmann constant in eV/K. return this._kB_evK = this._kB_evK || quantity('boltzmann', 'eV/K'); } } return new MemoizedConstants(); }36:[] 2e:{"title":"CosmoCalc","cells":[{"id":0,"value":"md`# CosmoCalc\n\nBasic calculations of background cosmology. To use:\n\\`\\`\\`\n import {Planck15,FlatLambdaCDM,LambdaCDM,w0waCDM} from '@dkirkby/cosmocalc'\n\\`\\`\\`\nThe class structure and API is similar to the python [astropy.cosmology](https://docs.astropy.org/en/stable/cosmology/) library.\n\nSee the [CosmoViz notebook](https://observablehq.com/@dkirkby/cosmoviz?collection=@dkirkby/cosmology) for some associated visualization utilities.\n\nTODO:\n - Handle \\`Om0=0\\` as a special case as ${tex`a\\rightarrow 0`} for the integrable singularity ${tex`\\propto a^{-1/2}`} in the horizon integral.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":1145,"value":"md`\nThe neutrino relative density is:\n${tex.block`\nf_\\nu(z) \\equiv \\frac{\\Omega_\\nu(z)}{\\Omega_\\gamma(z)} = \\frac{7}{8} \\left(\\frac{4}{11}\\right)^{4/3}\n\\frac{N_{eff}}{N_\\nu}\\, \\sum_{k=1}^{N_\\nu}\\, f\\left(\\frac{m_{\\nu,k}}{(1+z) k_B T_{\\nu,0}}\\right)\n`}\nwith ${tex`N_\\nu = \\text{floor}(N_{eff})`} and\n${tex.block`\nf(y) = \\frac{120}{7 \\pi^4} \\int_0^\\infty dx\\, \\frac{x^2 \\sqrt{x^2 + y^2}}{e^x + 1}\n`}\nWe approximate ${tex`f(y)`} using an analytical fitting formula, given in equation (26) of [Komatsu et al. 2011, ApJS 192, 18](https://iopscience.iop.org/article/10.1088/0067-0049/192/2/18):\n${tex.block`\nf(y) \\simeq \\left[ 1 + (A y)^p \\right]^{1/p}\n`}\nwith constants\n${tex.block`\nA = \\frac{180\\zeta(3)}{7 \\pi^4} \\simeq 0.3173 \\quad , \\quad p = 1.83 \\; .\n`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":982,"value":"md`Unrealistic models can have scales where ${tex`E(a)=0`} resulting in a local min/max of the expansion. These scales are solutions to the 4th-order polynomial:\n${tex.block`\n\\Omega_{\\Lambda,0} a^4 + \\left(1 - \\Omega_{r,0}-\\Omega_{m,0}-\\Omega_{\\Lambda,0}\\right) a^2\n+ \\Omega_{m,0} a + \\Omega_{r,0} = 0\n`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1147,"value":"md`\n*TODO:*\n - Add more derived quantity methods (angular/lumi distance, volume element, etc).\n - Calculate or input baryon fraction?\n - Implement better methods for selecting z or a integrand and integration order.\n - Cache efuncsq(z) and p=0,1 integrals for derived quantities to use.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":761,"value":"md`Calculate quantities for a Planck15 cosmology:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":755,"value":"Planck15.age // in Gyr","pinned":true,"hidden":false,"mode":"ojs"},{"id":1135,"value":"(1 - Planck15.lookback_time(0,4) / Planck15.age) * 50","pinned":true,"hidden":false,"mode":"ojs"},{"id":752,"value":"Planck15.critical_density0.toString()","pinned":true,"hidden":false,"mode":"ojs"},{"id":765,"value":"Planck15.Ode0","pinned":true,"hidden":false,"mode":"ojs"},{"id":940,"value":"Planck15.ngamma(0)","pinned":true,"hidden":false,"mode":"ojs"},{"id":966,"value":"Planck15.ndm(0,0.1) * 1e6","pinned":true,"hidden":false,"mode":"ojs"},{"id":950,"value":"Planck15.ndm(0) / Planck15.nb(0)","pinned":true,"hidden":false,"mode":"ojs"},{"id":866,"value":"{\n // z<0 corresponds to future times with negative comoving distances and lookback times.\n const z=linspace(-0.5, 2, 11, {dtype:Float64Array});\n return Planck15.tidyz(z, 'comoving_distance,lookback_time,critical_density,Ode');\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":833,"value":"Planck15._critical_density0","pinned":true,"hidden":false,"mode":"ojs"},{"id":692,"value":"md`Calculate quantities for a non-flat ${tex`\\Omega_{k,0}=0.1`} universe with ${tex`w_0=-0.9, w_a=0.1`}:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":615,"value":"w0wa = new w0waCDM({H0:70, Om0:0.25, Ok0:0.1, w0:-0.9, wa:0.1, Tcmb0:2.7, Neff:3.04})","pinned":true,"hidden":false,"mode":"ojs"},{"id":622,"value":"w0wa.age","pinned":true,"hidden":false,"mode":"ojs"},{"id":1113,"value":"md`Calculate quantities for pure fluid model with non-evolving dark-energy and neutrino sectors. Constructor args are ${tex`\\Omega_{r,0}, \\Omega_{m,0}, \\Omega_{de,0}, H_0`}.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1116,"value":"// Om0=0.3, Ode0=0.7 age in Gyr\nnew FluidsModel(0,0.3,0.7,70).age","pinned":true,"hidden":false,"mode":"ojs"},{"id":1127,"value":"// Model without radiation has infinite horizon\nnew FluidsModel(0,0.3,0.7,70).horizon","pinned":true,"hidden":false,"mode":"ojs"},{"id":1119,"value":"// Or0=Om0=1 universe reaches max size at a=1.6\nnew FluidsModel(1,1,0,70).aminmax","pinned":true,"hidden":false,"mode":"ojs"},{"id":1121,"value":"// Ode0=1 universe reaches min size at a=0.7\nnew FluidsModel(0,0,2,70).aminmax","pinned":true,"hidden":false,"mode":"ojs"},{"id":1252,"value":"Planck15NonAccel = Planck15.non_accelerated()","pinned":true,"hidden":false,"mode":"ojs"},{"id":701,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":163,"value":"Planck15 = new FlatLambdaCDM(Planck15Params);","pinned":false,"hidden":false,"mode":"ojs"},{"id":415,"value":"Planck15Params = ({H0:67.74, Om0:0.3075, Ob0:0.0486, Ok0:0, Tcmb0:2.7255, Neff:3.046, m_nu:[0, 0, 0.06], w0:-1, wa:0});","pinned":false,"hidden":false,"mode":"ojs"},{"id":382,"value":"class FlatLambdaCDM extends LambdaCDM {\n // Concrete subclass of LambdaCDM with zero curvature Ok0=0.\n // If Ok0 is present in the ctor params it must have the value Ok0.\n constructor(params={}) {\n if((params.Ok0 !== undefined) && (params.Ok0 != 0))\n throw Error('Cannot specify Ok0 != 0 for FlatLambdaCDM');\n super(Object.assign(params, {Ok0:0}));\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":362,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":396,"value":"class w0waCDM extends FLRW {\n constructor({H0,Om0,Ob0,Ok0,Tcmb0,Neff,m_nu,w0=-1,wa=0}={}) {\n super({H0:H0,Om0:Om0,Ob0:Ob0,Ok0:Ok0,Tcmb0:Tcmb0,Neff:Neff,m_nu:m_nu});\n this._w0 = w0;\n this._wa = wa;\n this._power = 3 * (1 + w0 + wa);\n this._coef = -3 * wa;\n }\n w(z) {\n return this._w0 + this._wa * z / (1 + z);\n }\n de_density_scale(z) {\n const zp1= 1 + z;\n return zp1**this._power * Math.exp(this._coef * z / zp1);\n }\n de_density_scale_a4(a) {\n // Calculate a**4 * de_density_scale(1/a-1).\n return a**(4-this._power) * Math.exp(this._coef * (1-a));\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":42,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":988,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":996,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":124,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":991,"value":"import {realsolns,quarticsolve} from '@dkirkby/poly'","pinned":false,"hidden":false,"mode":"ojs"},{"id":230,"value":"import {quad} from '@dkirkby/integrate'","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"import {quantity} from '@dkirkby/units'","pinned":false,"hidden":false,"mode":"ojs"},{"id":910,"value":"import {tidy,arraycat} from '@dkirkby/array-utilities'","pinned":false,"hidden":false,"mode":"ojs"},{"id":30,"value":"import {linspace,cumsum} from '@dkirkby/jumpy'","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W36","stdlib":"1","titling":"auto"}