import UIKit import SceneKit import PlaygroundSupport
// Constants let c: Double = 299792.458 // speed of light in km/s let H0: Double = 70 // Hubble constant in km/s/Mpc let G: Double = 6.6743e-11 // gravitational constant in m^3/kg/s^2 let rho_crit: Double = 1.878e-26 // critical density of the universe in kg/m^3
// Parameters let N: Int = 1000 // number of particles let t0: Double = 0 // initial time in Gyr let tf: Double = 13.8 // final time in Gyr let dt: Double = 0.01 // time step in Gyr let a0: Double = 1 / (1 + t0 * H0 / 2997.92) // initial scale factor
// Generate initial conditions let x = (0..<N).map { _ in return SCNVector3( CGFloat(Double.random(in: -10...10)), CGFloat(Double.random(in: -10...10)), CGFloat(Double.random(in: -10...10)) ) } let v = (0..<N).map { _ in return SCNVector3( CGFloat(Double.random(in: -100...100)), CGFloat(Double.random(in: -100...100)), CGFloat(Double.random(in: -100...100)) ) } let m = (0..<N).map { _ in return CGFloat(Double.random(in: 1e11...2e11)) }
// Define Friedmann equations func HubbleParameter(_ a: Double) -> Double { return H0 * sqrt((rho_crit / 3) * (1 / pow(a, 3))) }
func Acceleration(_ x: SCNVector3, _ v: SCNVector3) -> SCNVector3 { let r = x.length() let a = SCNVector3( -G * x.x / pow(r, 3), -G * x.y / pow(r, 3), -G * x.z / pow(r, 3) ) * m.reduce(0, +) a -= SCNVector3( Float(HubbleParameter(1) * Double(x.x)), Float(HubbleParameter(1) * Double(x.y)), Float(HubbleParameter(1) * Double(x.z)) ) return a }
// Initialize arrays for storing data let t = stride(from: t0, to: tf, by: dt).map { $0 } var a = Array(repeating: 0.0, count: t.count) var v = Array(repeating: SCNVector3Zero, count: t.count) var x = Array(repeating: SCNVector3Zero, count: t.count)
// Set initial conditions a[0] = a0 v[0] = v[0] x[0] = x[0]
// Integrate equations of motion using Verlet algorithm for i in 1..<t.count { x[i] = x[i-1] + v[i-1] * Float(dt) + Acceleration(x[i-1], v[i-1]) * pow(Float(dt), 2) / 2 a[i] = a0 * (1 + HubbleParameter(a0) * (t[i] - t0)) // scale factor as a function of time v[i] = v[i-1]
